MNE-CPP 0.1.9
A Framework for Electrophysiology
Loading...
Searching...
No Matches
mne_surface_or_volume.cpp
Go to the documentation of this file.
1//=============================================================================================================
37//=============================================================================================================
38// INCLUDES
39//=============================================================================================================
40
41#include "mne_surface_or_volume.h"
42#include "mne_surface_old.h"
43#include "mne_source_space_old.h"
44#include "mne_patch_info.h"
45//#include "fwd_bem_model.h"
46#include "mne_nearest.h"
47#include "filter_thread_arg.h"
48#include "mne_triangle.h"
49#include "mne_msh_display_surface.h"
50#include "mne_proj_data.h"
51#include "mne_vol_geom.h"
52#include "mne_mgh_tag_group.h"
53#include "mne_mgh_tag.h"
54
55#include <fiff/fiff_stream.h>
56#include <fiff/c/fiff_digitizer_data.h>
57#include <fiff/fiff_dig_point.h>
58
59#include <utils/sphere.h>
60#include <utils/ioutils.h>
61
62#include <QFile>
63#include <QCoreApplication>
64#include <QtConcurrent>
65
66#define _USE_MATH_DEFINES
67#include <math.h>
68
69#include <Eigen/Dense>
70#include <Eigen/Sparse>
71
72//ToDo don't use access and unlink -> use QT stuff instead
73#if defined(WIN32) || defined(_WIN32) || defined(__WIN32) && !defined(__CYGWIN__)
74#include <io.h>
75#else
76#include <unistd.h>
77#endif
78
79//=============================================================================================================
80// USED NAMESPACES
81//=============================================================================================================
82
83using namespace Eigen;
84using namespace FIFFLIB;
85using namespace MNELIB;
86
87//============================= dot.h =============================
88
89#ifndef TRUE
90#define TRUE 1
91#endif
92
93#ifndef FALSE
94#define FALSE 0
95#endif
96
97#ifndef FAIL
98#define FAIL -1
99#endif
100
101#ifndef OK
102#define OK 0
103#endif
104
105#define X_17 0
106#define Y_17 1
107#define Z_17 2
108
109#define VEC_DOT_17(x,y) ((x)[X_17]*(y)[X_17] + (x)[Y_17]*(y)[Y_17] + (x)[Z_17]*(y)[Z_17])
110
111#define VEC_LEN_17(x) sqrt(VEC_DOT_17(x,x))
112
113#define VEC_DIFF_17(from,to,diff) {\
114 (diff)[X_17] = (to)[X_17] - (from)[X_17];\
115 (diff)[Y_17] = (to)[Y_17] - (from)[Y_17];\
116 (diff)[Z_17] = (to)[Z_17] - (from)[Z_17];\
117 }
118
119#define VEC_COPY_17(to,from) {\
120 (to)[X_17] = (from)[X_17];\
121 (to)[Y_17] = (from)[Y_17];\
122 (to)[Z_17] = (from)[Z_17];\
123 }
124
125#define CROSS_PRODUCT_17(x,y,xy) {\
126 (xy)[X_17] = (x)[Y_17]*(y)[Z_17]-(y)[Y_17]*(x)[Z_17];\
127 (xy)[Y_17] = -((x)[X_17]*(y)[Z_17]-(y)[X_17]*(x)[Z_17]);\
128 (xy)[Z_17] = (x)[X_17]*(y)[Y_17]-(y)[X_17]*(x)[Y_17];\
129 }
130
131#define MALLOC_17(x,t) (t *)malloc((x)*sizeof(t))
132
133#define REALLOC_17(x,y,t) (t *)((x == NULL) ? malloc((y)*sizeof(t)) : realloc((x),(y)*sizeof(t)))
134
135#define ALLOC_INT_17(x) MALLOC_17(x,int)
136
137#define ALLOC_CMATRIX_17(x,y) mne_cmatrix_17((x),(y))
138
139static void matrix_error_17(int kind, int nr, int nc)
140
141{
142 if (kind == 1)
143 printf("Failed to allocate memory pointers for a %d x %d matrix\n",nr,nc);
144 else if (kind == 2)
145 printf("Failed to allocate memory for a %d x %d matrix\n",nr,nc);
146 else
147 printf("Allocation error for a %d x %d matrix\n",nr,nc);
148 if (sizeof(void *) == 4) {
149 printf("This is probably because you seem to be using a computer with 32-bit architecture.\n");
150 printf("Please consider moving to a 64-bit platform.");
151 }
152 printf("Cannot continue. Sorry.\n");
153 exit(1);
154}
155
156float** mne_cmatrix_17(int numPoints,int numDim)
157{
158 float** m;
159 float* whole;
160
161 m = MALLOC_17(numPoints, float *);
162 if (!m)
163 {
164 matrix_error_17( 1, numPoints, numDim);
165 }
166
167 whole = MALLOC_17( numPoints * numDim, float);
168 if (!whole)
169 {
170 matrix_error_17(2, numPoints, numDim);
171 }
172
173 for(int i = 0; i < numPoints; ++i)
174 {
175 m[i] = &whole[ i * numDim ];
176 }
177
178 return m;
179}
180
181#define FREE_17(x) if ((char *)(x) != NULL) free((char *)(x))
182
183#define FREE_CMATRIX_17(m) mne_free_cmatrix_17((m))
184
185#define FREE_ICMATRIX_17(m) mne_free_icmatrix_17((m))
186
187void mne_free_cmatrix_17 (float **m)
188{
189 if (m) {
190 FREE_17(*m);
191 FREE_17(m);
192 }
193}
194
195void mne_free_icmatrix_17 (int **m)
196
197{
198 if (m) {
199 FREE_17(*m);
200 FREE_17(m);
201 }
202}
203
204#define NNEIGHBORS 26
205
206#define CURVATURE_FILE_MAGIC_NUMBER (16777215)
207
208#define TAG_MGH_XFORM 31
209#define TAG_SURF_GEOM 21
210#define TAG_OLD_USEREALRAS 2
211#define TAG_COLORTABLE 5
212#define TAG_OLD_MGH_XFORM 30
213#define TAG_OLD_COLORTABLE 1
214
215#define TAG_USEREALRAS 4
216
217#define ALLOC_ICMATRIX_17(x,y) mne_imatrix_17((x),(y))
218
219int **mne_imatrix_17(int nr,int nc)
220
221{
222 int i,**m;
223 int *whole;
224
225 m = MALLOC_17(nr,int *);
226 if (!m) matrix_error_17(1,nr,nc);
227 whole = MALLOC_17(nr*nc,int);
228 if (!whole) matrix_error_17(2,nr,nc);
229 for(i=0;i<nr;i++)
230 m[i] = whole + i*nc;
231 return m;
232}
233
234//float
235Eigen::MatrixXf toFloatEigenMatrix_17(float **mat, const int m, const int n)
236{
237 Eigen::MatrixXf eigen_mat(m,n);
238
239 for ( int i = 0; i < m; ++i)
240 for ( int j = 0; j < n; ++j)
241 eigen_mat(i,j) = mat[i][j];
242
243 return eigen_mat;
244}
245
246void fromFloatEigenMatrix_17(const Eigen::MatrixXf& from_mat, float **& to_mat, const int m, const int n)
247{
248 for ( int i = 0; i < m; ++i)
249 for ( int j = 0; j < n; ++j)
250 to_mat[i][j] = from_mat(i,j);
251}
252
253void fromFloatEigenMatrix_17(const Eigen::MatrixXf& from_mat, float **& to_mat)
254{
255 fromFloatEigenMatrix_17(from_mat, to_mat, from_mat.rows(), from_mat.cols());
256}
257
258void fromIntEigenMatrix_17(const Eigen::MatrixXi& from_mat, int **&to_mat, const int m, const int n)
259{
260 for ( int i = 0; i < m; ++i)
261 for ( int j = 0; j < n; ++j)
262 to_mat[i][j] = from_mat(i,j);
263}
264
265void fromIntEigenMatrix_17(const Eigen::MatrixXi& from_mat, int **&to_mat)
266{
267 fromIntEigenMatrix_17(from_mat, to_mat, from_mat.rows(), from_mat.cols());
268}
269
270//============================= make_volume_source_space.c =============================
271
272static FiffCoordTransOld* make_voxel_ras_trans(float *r0,
273 float *x_ras,
274 float *y_ras,
275 float *z_ras,
276 float *voxel_size)
277
278{
279 FiffCoordTransOld* t;
280 float rot[3][3],move[3];
281 int j,k;
282
283 VEC_COPY_17(move,r0);
284
285 for (j = 0; j < 3; j++) {
286 rot[j][0] = x_ras[j];
287 rot[j][1] = y_ras[j];
288 rot[j][2] = z_ras[j];
289 }
290
291 for (j = 0; j < 3; j++)
292 for (k = 0; k < 3; k++)
293 rot[j][k] = voxel_size[k]*rot[j][k];
294
295 t = new FiffCoordTransOld(FIFFV_MNE_COORD_MRI_VOXEL,FIFFV_COORD_MRI,rot,move);
296
297 return t;
298}
299
300static int comp_points1(const void *vp1,const void *vp2)
301
302{
303 MneNearest* v1 = (MneNearest*)vp1;
304 MneNearest* v2 = (MneNearest*)vp2;
305
306 if (v1->nearest > v2->nearest)
307 return 1;
308 else if (v1->nearest == v2->nearest)
309 return 0;
310 else
311 return -1;
312}
313
314static int comp_points2(const void *vp1,const void *vp2)
315
316{
317 MneNearest* v1 = (MneNearest*)vp1;
318 MneNearest* v2 = (MneNearest*)vp2;
319
320 if (v1->vert > v2->vert)
321 return 1;
322 else if (v1->vert == v2->vert)
323 return 0;
324 else
325 return -1;
326}
327
328void mne_sort_nearest_by_nearest(MneNearest* points, int npoint)
329
330{
331 if (npoint > 1 && points != NULL)
332 qsort(points,npoint,sizeof(MneNearest),comp_points1);
333 return;
334}
335
336//=============================================================================================================
337// DEFINE MEMBER METHODS
338//=============================================================================================================
339
343
344//=============================================================================================================
345
346MneSurfaceOrVolume::~MneSurfaceOrVolume()
347{
348 int k;
349 FREE_CMATRIX_17(this->rr);
350 FREE_CMATRIX_17(this->nn);
351 FREE_17(this->inuse);
352 FREE_17(this->vertno);
353 FREE_17(this->tris);
354 FREE_ICMATRIX_17(this->itris);
355
356 FREE_17(this->use_tris);
357 FREE_ICMATRIX_17(this->use_itris);
358 if (this->neighbor_tri) {
359 for (k = 0; k < this->np; k++)
360 FREE_17(this->neighbor_tri[k]);
361 FREE_17(this->neighbor_tri);
362 }
363 FREE_17(this->nneighbor_tri);
364 FREE_17(this->curv);
365
366 if (this->neighbor_vert) {
367 for (k = 0; k < this->np; k++)
368 FREE_17(this->neighbor_vert[k]);
369 FREE_17(this->neighbor_vert);
370 }
371 FREE_17(this->nneighbor_vert);
372 if (this->vert_dist) {
373 for (k = 0; k < this->np; k++)
374 FREE_17(this->vert_dist[k]);
375 FREE_17(this->vert_dist);
376 }
377 FREE_17(this->nearest);
378 if (this->patches) {
379 for (k = 0; k < this->npatch; k++)
380 if(this->patches[k])
381 delete this->patches[k];
382 FREE_17(this->patches);
383 }
384 if(this->dist)
385 delete this->dist;
386 FREE_17(this->voxel_surf_RAS_t);
387 FREE_17(this->MRI_voxel_surf_RAS_t);
388 FREE_17(this->MRI_surf_RAS_RAS_t);
389 if(this->interpolator)
390 delete this->interpolator;
391 this->MRI_volume.clear();
392
393 if(this->vol_geom)
394 delete this->vol_geom;
395 delete((MneMghTagGroup*)this->mgh_tags);
396
397 if (this->user_data && this->user_data_free)
398 this->user_data_free(this->user_data);
399}
400
401//=============================================================================================================
402
403double MneSurfaceOrVolume::solid_angle(float *from, MneTriangle* tri) /* ...to this triangle */
404/*
405 * Compute the solid angle according to van Oosterom's
406 * formula
407 */
408{
409 double v1[3],v2[3],v3[3];
410 double l1,l2,l3,s,triple;
411 double cross[3];
412
413 VEC_DIFF_17 (from,tri->r1,v1);
414 VEC_DIFF_17 (from,tri->r2,v2);
415 VEC_DIFF_17 (from,tri->r3,v3);
416
417 CROSS_PRODUCT_17(v1,v2,cross);
418 triple = VEC_DOT_17(cross,v3);
419
420 l1 = VEC_LEN_17(v1);
421 l2 = VEC_LEN_17(v2);
422 l3 = VEC_LEN_17(v3);
423 s = (l1*l2*l3+VEC_DOT_17(v1,v2)*l3+VEC_DOT_17(v1,v3)*l2+VEC_DOT_17(v2,v3)*l1);
424
425 return (2.0*atan2(triple,s));
426}
427
428//=============================================================================================================
429
430double MneSurfaceOrVolume::sum_solids(float *from, MneSurfaceOld* surf)
431{
432 int k;
433 double tot_angle, angle;
434 for (k = 0, tot_angle = 0.0; k < surf->ntri; k++) {
435 angle = solid_angle(from,surf->tris+k);
436 tot_angle += angle;
437 }
438 return tot_angle;
439}
440
441//=============================================================================================================
442
443int MneSurfaceOrVolume::mne_filter_source_spaces(MneSurfaceOld* surf, float limit, FiffCoordTransOld* mri_head_t, MneSourceSpaceOld* *spaces, int nspace, FILE *filtered) /* Provide a list of filtered points here */
444/*
445 * Remove all source space points closer to the surface than a given limit
446 */
447{
448 MneSourceSpaceOld* s;
449 int k,p1,p2;
450 float r1[3];
451 float mindist,dist,diff[3];
452 int minnode;
453 int omit,omit_outside;
454 double tot_angle;
455
456 if (surf == NULL)
457 return OK;
458 if (spaces[0]->coord_frame == FIFFV_COORD_HEAD && mri_head_t == NULL) {
459 printf("Source spaces are in head coordinates and no coordinate transform was provided!");
460 return FAIL;
461 }
462 /*
463 * How close are the source points to the surface?
464 */
465 printf("Source spaces are in ");
466 if (spaces[0]->coord_frame == FIFFV_COORD_HEAD)
467 printf("head coordinates.\n");
468 else if (spaces[0]->coord_frame == FIFFV_COORD_MRI)
469 printf("MRI coordinates.\n");
470 else
471 printf("unknown (%d) coordinates.\n",spaces[0]->coord_frame);
472 printf("Checking that the sources are inside the bounding surface ");
473 if (limit > 0.0)
474 printf("and at least %6.1f mm away",1000*limit);
475 printf(" (will take a few...)\n");
476 omit = 0;
477 omit_outside = 0;
478 for (k = 0; k < nspace; k++) {
479 s = spaces[k];
480 for (p1 = 0; p1 < s->np; p1++)
481 if (s->inuse[p1]) {
482 VEC_COPY_17(r1,s->rr[p1]); /* Transform the point to MRI coordinates */
483 if (s->coord_frame == FIFFV_COORD_HEAD)
484 FiffCoordTransOld::fiff_coord_trans_inv(r1,mri_head_t,FIFFV_MOVE);
485 /*
486 * Check that the source is inside the inner skull surface
487 */
488 tot_angle = sum_solids(r1,surf)/(4*M_PI);
489 if (std::fabs(tot_angle-1.0) > 1e-5) {
490 omit_outside++;
491 s->inuse[p1] = FALSE;
492 s->nuse--;
493 if (filtered)
494 fprintf(filtered,"%10.3f %10.3f %10.3f\n",
495 1000*r1[X_17],1000*r1[Y_17],1000*r1[Z_17]);
496 }
497 else if (limit > 0.0) {
498 /*
499 * Check the distance limit
500 */
501 mindist = 1.0;
502 minnode = 0;
503 for (p2 = 0; p2 < surf->np; p2++) {
504 VEC_DIFF_17(r1,surf->rr[p2],diff);
505 dist = VEC_LEN_17(diff);
506 if (dist < mindist) {
507 mindist = dist;
508 minnode = p2;
509 }
510 }
511 if (mindist < limit) {
512 omit++;
513 s->inuse[p1] = FALSE;
514 s->nuse--;
515 if (filtered)
516 fprintf(filtered,"%10.3f %10.3f %10.3f\n",
517 1000*r1[X_17],1000*r1[Y_17],1000*r1[Z_17]);
518 }
519 }
520 }
521 }
522 (void)minnode; // squash compiler warning, this is unused
523 if (omit_outside > 0)
524 printf("%d source space points omitted because they are outside the inner skull surface.\n",
525 omit_outside);
526 if (omit > 0)
527 printf("%d source space points omitted because of the %6.1f-mm distance limit.\n",
528 omit,1000*limit);
529 printf("Thank you for waiting.\n");
530 return OK;
531}
532
533//=============================================================================================================
534
535int MneSurfaceOrVolume::mne_add_patch_stats(MneSourceSpaceOld* s)
536{
537 MneNearest* nearest = s->nearest;
538 MneNearest* this_patch;
539 MnePatchInfo* *pinfo = MALLOC_17(s->nuse,MnePatchInfo*);
540 int nave,p,q,k;
541
542 printf("Computing patch statistics...\n");
543 if (!s->neighbor_tri)
544 if (mne_source_space_add_geometry_info(s,FALSE) != OK)
545 goto bad;
546
547 if (s->nearest == NULL) {
548 qCritical("The patch information is not available.");
549 goto bad;
550 }
551 if (s->nuse == 0) {
552 FREE_17(s->patches);
553 s->patches = NULL;
554 s->npatch = 0;
555 return OK;
556 }
557 /*
558 * Calculate the average normals and the patch areas
559 */
560 printf("\tareas, average normals, and mean deviations...");
561 mne_sort_nearest_by_nearest(nearest,s->np);
562 nave = 1;
563 for (p = 1, q = 0; p < s->np; p++) {
564 if (nearest[p].nearest != nearest[p-1].nearest) {
565 if (nave == 0) {
566 qCritical("No vertices belong to the patch of vertex %d",nearest[p-1].nearest);
567 goto bad;
568 }
569 if (s->vertno[q] == nearest[p-1].nearest) { /* Some source space points may have been omitted since
570 * the patch information was computed */
571 pinfo[q] = new MnePatchInfo();
572 pinfo[q]->vert = nearest[p-1].nearest;
573 this_patch = nearest+p-nave;
574 pinfo[q]->memb_vert = MALLOC_17(nave,int);
575 pinfo[q]->nmemb = nave;
576 for (k = 0; k < nave; k++) {
577 pinfo[q]->memb_vert[k] = this_patch[k].vert;
578 this_patch[k].patch = pinfo[q];
579 }
580 MnePatchInfo::calculate_patch_area(s,pinfo[q]);
581 MnePatchInfo::calculate_normal_stats(s,pinfo[q]);
582 q++;
583 }
584 nave = 0;
585 }
586 nave++;
587 }
588 if (nave == 0) {
589 qCritical("No vertices belong to the patch of vertex %d",nearest[p-1].nearest);
590 goto bad;
591 }
592 if (s->vertno[q] == nearest[p-1].nearest) {
593 pinfo[q] = new MnePatchInfo;
594 pinfo[q]->vert = nearest[p-1].nearest;
595 this_patch = nearest+p-nave;
596 pinfo[q]->memb_vert = MALLOC_17(nave,int);
597 pinfo[q]->nmemb = nave;
598 for (k = 0; k < nave; k++) {
599 pinfo[q]->memb_vert[k] = this_patch[k].vert;
600 this_patch[k].patch = pinfo[q];
601 }
602 MnePatchInfo::calculate_patch_area(s,pinfo[q]);
603 MnePatchInfo::calculate_normal_stats(s,pinfo[q]);
604 q++;
605 }
606 printf(" %d/%d [done]\n",q,s->nuse);
607
608 if (s->patches) {
609 for (k = 0; k < s->npatch; k++)
610 if(s->patches[k])
611 delete s->patches[k];
612 FREE_17(s->patches);
613 }
614 s->patches = pinfo;
615 s->npatch = s->nuse;
616
617 return OK;
618
619bad : {
620 FREE_17(pinfo);
621 return FAIL;
622 }
623}
624
625//=============================================================================================================
626
627void MneSurfaceOrVolume::rearrange_source_space(MneSourceSpaceOld* s)
628{
629 int k,p;
630
631 for (k = 0, s->nuse = 0; k < s->np; k++)
632 if (s->inuse[k])
633 s->nuse++;
634
635 if (s->nuse == 0) {
636 FREE_17(s->vertno);
637 s->vertno = NULL;
638 }
639 else {
640 s->vertno = REALLOC_17(s->vertno,s->nuse,int);
641 for (k = 0, p = 0; k < s->np; k++)
642 if (s->inuse[k])
643 s->vertno[p++] = k;
644 }
645 if (s->nearest)
646 mne_add_patch_stats(s);
647 return;
648}
649
650//=============================================================================================================
651
652void *MneSurfaceOrVolume::filter_source_space(void *arg)
653{
654 FilterThreadArg* a = (FilterThreadArg*)arg;
655 int p1,p2;
656 double tot_angle;
657 int omit,omit_outside;
658 float r1[3];
659 float mindist,dist,diff[3];
660 int minnode;
661
662 omit = 0;
663 omit_outside = 0;
664
665 for (p1 = 0; p1 < a->s->np; p1++) {
666 if (a->s->inuse[p1]) {
667 VEC_COPY_17(r1,a->s->rr[p1]); /* Transform the point to MRI coordinates */
668 if (a->s->coord_frame == FIFFV_COORD_HEAD)
669 FiffCoordTransOld::fiff_coord_trans_inv(r1,a->mri_head_t,FIFFV_MOVE);
670 /*
671 * Check that the source is inside the inner skull surface
672 */
673 tot_angle = sum_solids(r1,a->surf)/(4*M_PI);
674 if (std::fabs(tot_angle-1.0) > 1e-5) {
675 omit_outside++;
676 a->s->inuse[p1] = FALSE;
677 a->s->nuse--;
678 if (a->filtered)
679 fprintf(a->filtered,"%10.3f %10.3f %10.3f\n",
680 1000*r1[X_17],1000*r1[Y_17],1000*r1[Z_17]);
681 }
682 else if (a->limit > 0.0) {
683 /*
684 * Check the distance limit
685 */
686 mindist = 1.0;
687 minnode = 0;
688 for (p2 = 0; p2 < a->surf->np; p2++) {
689 VEC_DIFF_17(r1,a->surf->rr[p2],diff);
690 dist = VEC_LEN_17(diff);
691 if (dist < mindist) {
692 mindist = dist;
693 minnode = p2;
694 }
695 }
696 if (mindist < a->limit) {
697 omit++;
698 a->s->inuse[p1] = FALSE;
699 a->s->nuse--;
700 if (a->filtered)
701 fprintf(a->filtered,"%10.3f %10.3f %10.3f\n",
702 1000*r1[X_17],1000*r1[Y_17],1000*r1[Z_17]);
703 }
704 }
705 }
706 }
707 (void)minnode; // squash compiler warning, set but unused
708 if (omit_outside > 0)
709 printf("%d source space points omitted because they are outside the inner skull surface.\n",
710 omit_outside);
711 if (omit > 0)
712 printf("%d source space points omitted because of the %6.1f-mm distance limit.\n",
713 omit,1000*a->limit);
714 a->stat = OK;
715 return NULL;
716}
717
718//=============================================================================================================
719
720int MneSurfaceOrVolume::filter_source_spaces(float limit, char *bemfile, FiffCoordTransOld *mri_head_t, MneSourceSpaceOld* *spaces, int nspace, FILE *filtered, bool use_threads) /* Use multiple threads if possible? */
721/*
722 * Remove all source space points closer to the surface than a given limit
723 */
724{
725 MneSurfaceOld* surf = NULL;
726 int k;
727 int nproc = QThread::idealThreadCount();
728 FilterThreadArg* a;
729
730 if (!bemfile)
731 return OK;
732
733 if ((surf = MneSurfaceOld::read_bem_surface(bemfile,FIFFV_BEM_SURF_ID_BRAIN,FALSE,NULL)) == NULL) {
734 qCritical("BEM model does not have the inner skull triangulation!");
735 return FAIL;
736 }
737 /*
738 * How close are the source points to the surface?
739 */
740 printf("Source spaces are in ");
741 if (spaces[0]->coord_frame == FIFFV_COORD_HEAD)
742 printf("head coordinates.\n");
743 else if (spaces[0]->coord_frame == FIFFV_COORD_MRI)
744 printf("MRI coordinates.\n");
745 else
746 printf("unknown (%d) coordinates.\n",spaces[0]->coord_frame);
747 printf("Checking that the sources are inside the inner skull ");
748 if (limit > 0.0)
749 printf("and at least %6.1f mm away",1000*limit);
750 printf(" (will take a few...)\n");
751 if (nproc < 2 || nspace == 1 || !use_threads) {
752 /*
753 * This is the conventional calculation
754 */
755 for (k = 0; k < nspace; k++) {
756 a = new FilterThreadArg();
757 a->s = spaces[k];
758 a->mri_head_t = mri_head_t;
759 a->surf = surf;
760 a->limit = limit;
761 a->filtered = filtered;
762 filter_source_space(a);
763 if(a)
764 delete a;
765 rearrange_source_space(spaces[k]);
766 }
767 }
768 else {
769 /*
770 * Calculate all (both) source spaces simultaneously
771 */
772 QList<FilterThreadArg*> args;//filterThreadArg *args = MALLOC_17(nspace,filterThreadArg);
773
774 for (k = 0; k < nspace; k++) {
775 a = new FilterThreadArg();
776 a->s = spaces[k];
777 a->mri_head_t = mri_head_t;
778 a->surf = surf;
779 a->limit = limit;
780 a->filtered = filtered;
781 args.append(a);
782 }
783 /*
784 * Ready to start the threads & Wait for them to complete
785 */
786 QtConcurrent::blockingMap(args, filter_source_space);
787
788 for (k = 0; k < nspace; k++) {
789 rearrange_source_space(spaces[k]);
790 if(args[k])
791 delete args[k];
792 }
793 }
794 if(surf)
795 delete surf;
796 printf("Thank you for waiting.\n\n");
797
798 return OK;
799}
800
801//=============================================================================================================
802
803MneSourceSpaceOld* MneSurfaceOrVolume::make_volume_source_space(MneSurfaceOld* surf, float grid, float exclude, float mindist)
804/*
805 * Make a source space which covers the volume bounded by surf
806 */
807{
808 float min[3],max[3],cm[3];
809 int minn[3],maxn[3];
810 float *node,maxdist,dist,diff[3];
811 int k,c;
812 MneSourceSpaceOld* sp = NULL;
813 int np,nplane,nrow;
814 int *neigh,nneigh;
815 int x,y,z;
816 /*
817 * Figure out the grid size
818 */
819 cm[X_17] = cm[Y_17] = cm[Z_17] = 0.0;
820 node = surf->rr[0];
821 for (c = 0; c < 3; c++)
822 min[c] = max[c] = node[c];
823
824 for (k = 0; k < surf->np; k++) {
825 node = surf->rr[k];
826 for (c = 0; c < 3; c++) {
827 cm[c] += node[c];
828 if (node[c] < min[c])
829 min[c] = node[c];
830 if (node[c] > max[c])
831 max[c] = node[c];
832 }
833 }
834 for (c = 0; c < 3; c++)
835 cm[c] = cm[c]/surf->np;
836 /*
837 * Define the sphere which fits the surface
838 */
839 maxdist = 0.0;
840 for (k = 0; k < surf->np; k++) {
841 VEC_DIFF_17(cm,surf->rr[k],diff);
842 dist = VEC_LEN_17(diff);
843 if (dist > maxdist)
844 maxdist = dist;
845 }
846 printf("Surface CM = (%6.1f %6.1f %6.1f) mm\n",
847 1000*cm[X_17], 1000*cm[Y_17], 1000*cm[Z_17]);
848 printf("Surface fits inside a sphere with radius %6.1f mm\n",1000*maxdist);
849 printf("Surface extent:\n"
850 "\tx = %6.1f ... %6.1f mm\n"
851 "\ty = %6.1f ... %6.1f mm\n"
852 "\tz = %6.1f ... %6.1f mm\n",
853 1000*min[X_17],1000*max[X_17],
854 1000*min[Y_17],1000*max[Y_17],
855 1000*min[Z_17],1000*max[Z_17]);
856 for (c = 0; c < 3; c++) {
857 if (max[c] > 0)
858 maxn[c] = floor(std::fabs(max[c])/grid)+1;
859 else
860 maxn[c] = -floor(std::fabs(max[c])/grid)-1;
861 if (min[c] > 0)
862 minn[c] = floor(std::fabs(min[c])/grid)+1;
863 else
864 minn[c] = -floor(std::fabs(min[c])/grid)-1;
865 }
866 printf("Grid extent:\n"
867 "\tx = %6.1f ... %6.1f mm\n"
868 "\ty = %6.1f ... %6.1f mm\n"
869 "\tz = %6.1f ... %6.1f mm\n",
870 1000*(minn[X_17]*grid),1000*(maxn[X_17]*grid),
871 1000*(minn[Y_17]*grid),1000*(maxn[Y_17]*grid),
872 1000*(minn[Z_17]*grid),1000*(maxn[Z_17]*grid));
873 /*
874 * Now make the initial grid
875 */
876 np = 1;
877 for (c = 0; c < 3; c++)
878 np = np*(maxn[c]-minn[c]+1);
879 nplane = (maxn[X_17]-minn[X_17]+1)*(maxn[Y_17]-minn[Y_17]+1);
880 nrow = (maxn[X_17]-minn[X_17]+1);
881 sp = MneSurfaceOrVolume::mne_new_source_space(np);
882 sp->type = MNE_SOURCE_SPACE_VOLUME;
883 sp->nneighbor_vert = MALLOC_17(sp->np,int);
884 sp->neighbor_vert = MALLOC_17(sp->np,int *);
885 for (k = 0; k < sp->np; k++) {
886 sp->inuse[k] = TRUE;
887 sp->vertno[k] = k;
888 sp->nn[k][X_17] = sp->nn[k][Y_17] = 0.0; /* Source orientation is immaterial */
889 sp->nn[k][Z_17] = 1.0;
890 sp->neighbor_vert[k] = neigh = MALLOC_17(NNEIGHBORS,int);
891 sp->nneighbor_vert[k] = nneigh = NNEIGHBORS;
892 for (c = 0; c < nneigh; c++)
893 neigh[c] = -1;
894 sp->nuse++;
895 }
896 for (k = 0, z = minn[Z_17]; z <= maxn[Z_17]; z++) {
897 for (y = minn[Y_17]; y <= maxn[Y_17]; y++) {
898 for (x = minn[X_17]; x <= maxn[X_17]; x++, k++) {
899 sp->rr[k][X_17] = x*grid;
900 sp->rr[k][Y_17] = y*grid;
901 sp->rr[k][Z_17] = z*grid;
902 /*
903 * Figure out the neighborhood:
904 * 6-neighborhood first
905 */
906 neigh = sp->neighbor_vert[k];
907 if (z > minn[Z_17])
908 neigh[0] = k - nplane;
909 if (x < maxn[X_17])
910 neigh[1] = k + 1;
911 if (y < maxn[Y_17])
912 neigh[2] = k + nrow;
913 if (x > minn[X_17])
914 neigh[3] = k - 1;
915 if (y > minn[Y_17])
916 neigh[4] = k - nrow;
917 if (z < maxn[Z_17])
918 neigh[5] = k + nplane;
919 /*
920 * Then the rest to complete the 26-neighborhood
921 * First the plane below
922 */
923 if (z > minn[Z_17]) {
924 if (x < maxn[X_17]) {
925 neigh[6] = k + 1 - nplane;
926 if (y < maxn[Y_17])
927 neigh[7] = k + 1 + nrow - nplane;
928 }
929 if (y < maxn[Y_17])
930 neigh[8] = k + nrow - nplane;
931 if (x > minn[X_17]) {
932 if (y < maxn[Y_17])
933 neigh[9] = k - 1 + nrow - nplane;
934 neigh[10] = k - 1 - nplane;
935 if (y > minn[Y_17])
936 neigh[11] = k - 1 - nrow - nplane;
937 }
938 if (y > minn[Y_17]) {
939 neigh[12] = k - nrow - nplane;
940 if (x < maxn[X_17])
941 neigh[13] = k + 1 - nrow - nplane;
942 }
943 }
944 /*
945 * Then the same plane
946 */
947 if (x < maxn[X_17] && y < maxn[Y_17])
948 neigh[14] = k + 1 + nrow;
949 if (x > minn[X_17]) {
950 if (y < maxn[Y_17])
951 neigh[15] = k - 1 + nrow;
952 if (y > minn[Y_17])
953 neigh[16] = k - 1 - nrow;
954 }
955 if (y > minn[Y_17] && x < maxn[X_17])
956 neigh[17] = k + 1 - nrow - nplane;
957 /*
958 * Finally one plane above
959 */
960 if (z < maxn[Z_17]) {
961 if (x < maxn[X_17]) {
962 neigh[18] = k + 1 + nplane;
963 if (y < maxn[Y_17])
964 neigh[19] = k + 1 + nrow + nplane;
965 }
966 if (y < maxn[Y_17])
967 neigh[20] = k + nrow + nplane;
968 if (x > minn[X_17]) {
969 if (y < maxn[Y_17])
970 neigh[21] = k - 1 + nrow + nplane;
971 neigh[22] = k - 1 + nplane;
972 if (y > minn[Y_17])
973 neigh[23] = k - 1 - nrow + nplane;
974 }
975 if (y > minn[Y_17]) {
976 neigh[24] = k - nrow + nplane;
977 if (x < maxn[X_17])
978 neigh[25] = k + 1 - nrow + nplane;
979 }
980 }
981 }
982 }
983 }
984 printf("%d sources before omitting any.\n",sp->nuse);
985 /*
986 * Exclude infeasible points
987 */
988 for (k = 0; k < sp->np; k++) {
989 VEC_DIFF_17(cm,sp->rr[k],diff);
990 dist = VEC_LEN_17(diff);
991 if (dist < exclude || dist > maxdist) {
992 sp->inuse[k] = FALSE;
993 sp->nuse--;
994 }
995 }
996 printf("%d sources after omitting infeasible sources.\n",sp->nuse);
997 if (mne_filter_source_spaces(surf,mindist,NULL,&sp,1,NULL) != OK)
998 goto bad;
999 printf("%d sources remaining after excluding the sources outside the surface and less than %6.1f mm inside.\n",sp->nuse,1000*mindist);
1000 /*
1001 * Omit unused vertices from the neighborhoods
1002 */
1003 printf("Adjusting the neighborhood info...");
1004 for (k = 0; k < sp->np; k++) {
1005 neigh = sp->neighbor_vert[k];
1006 nneigh = sp->nneighbor_vert[k];
1007 if (sp->inuse[k]) {
1008 for (c = 0; c < nneigh; c++)
1009 if (!sp->inuse[neigh[c]])
1010 neigh[c] = -1;
1011 }
1012 else {
1013 for (c = 0; c < nneigh; c++)
1014 neigh[c] = -1;
1015 }
1016 }
1017 printf("[done]\n");
1018 /*
1019 * Set up the volume data (needed for creating the interpolation matrix)
1020 */
1021 {
1022 float r0[3],voxel_size[3],x_ras[3],y_ras[3],z_ras[3];
1023 int width,height,depth;
1024
1025 r0[X_17] = minn[X_17]*grid;
1026 r0[Y_17] = minn[Y_17]*grid;
1027 r0[Z_17] = minn[Z_17]*grid;
1028
1029 voxel_size[0] = grid;
1030 voxel_size[1] = grid;
1031 voxel_size[2] = grid;
1032
1033 width = (maxn[X_17]-minn[X_17]+1);
1034 height = (maxn[Y_17]-minn[Y_17]+1);
1035 depth = (maxn[Z_17]-minn[Z_17]+1);
1036
1037 for (k = 0; k < 3; k++)
1038 x_ras[k] = y_ras[k] = z_ras[k] = 0.0;
1039
1040 x_ras[0] = 1.0;
1041 y_ras[1] = 1.0;
1042 z_ras[2] = 1.0;
1043
1044 if ((sp->voxel_surf_RAS_t = make_voxel_ras_trans(r0,x_ras,y_ras,z_ras,voxel_size)) == NULL)
1045 goto bad;
1046
1047 sp->vol_dims[0] = width;
1048 sp->vol_dims[1] = height;
1049 sp->vol_dims[2] = depth;
1050 VEC_COPY_17(sp->voxel_size,voxel_size);
1051 }
1052
1053 return sp;
1054
1055bad : {
1056 if(sp)
1057 delete sp;
1058 return NULL;
1059 }
1060}
1061
1062//=============================================================================================================
1063
1064MneSourceSpaceOld* MneSurfaceOrVolume::mne_new_source_space(int np)
1065/*
1066 * Create a new source space and all associated data
1067 */
1068{
1069 MneSourceSpaceOld* res = new MneSourceSpaceOld();
1070 res->np = np;
1071 if (np > 0) {
1072 res->rr = ALLOC_CMATRIX_17(np,3);
1073 res->nn = ALLOC_CMATRIX_17(np,3);
1074 res->inuse = ALLOC_INT_17(np);
1075 res->vertno = ALLOC_INT_17(np);
1076 }
1077 else {
1078 res->rr = NULL;
1079 res->nn = NULL;
1080 res->inuse = NULL;
1081 res->vertno = NULL;
1082 }
1083 res->nuse = 0;
1084 res->ntri = 0;
1085 res->tris = NULL;
1086 res->itris = NULL;
1087 res->tot_area = 0.0;
1088
1089 res->nuse_tri = 0;
1090 res->use_tris = NULL;
1091 res->use_itris = NULL;
1092
1093 res->neighbor_tri = NULL;
1094 res->nneighbor_tri = NULL;
1095 res->curv = NULL;
1096 res->val = NULL;
1097
1098 res->neighbor_vert = NULL;
1099 res->nneighbor_vert = NULL;
1100 res->vert_dist = NULL;
1101
1102 res->coord_frame = FIFFV_COORD_MRI;
1103 res->id = FIFFV_MNE_SURF_UNKNOWN;
1104 res->subject.clear();
1105 res->type = FIFFV_MNE_SPACE_SURFACE;
1106
1107 res->nearest = NULL;
1108 res->patches = NULL;
1109 res->npatch = 0;
1110
1111 res->dist = NULL;
1112 res->dist_limit = -1.0;
1113
1114 res->voxel_surf_RAS_t = NULL;
1115 res->vol_dims[0] = res->vol_dims[1] = res->vol_dims[2] = 0;
1116
1117 res->MRI_volume.clear();
1118 res->MRI_surf_RAS_RAS_t = NULL;
1119 res->MRI_voxel_surf_RAS_t = NULL;
1120 res->MRI_vol_dims[0] = res->MRI_vol_dims[1] = res->MRI_vol_dims[2] = 0;
1121 res->interpolator = NULL;
1122
1123 res->vol_geom = NULL;
1124 res->mgh_tags = NULL;
1125 res->user_data = NULL;
1126 res->user_data_free = NULL;
1127
1128 res->cm[X_17] = res->cm[Y_17] = res->cm[Z_17] = 0.0;
1129
1130 return res;
1131}
1132
1133//=============================================================================================================
1134
1135MneSurfaceOld* MneSurfaceOrVolume::read_bem_surface(const QString &name, int which, int add_geometry, float *sigmap) /* Conductivity? */
1136{
1137 return read_bem_surface(name,which,add_geometry,sigmap,true);
1138}
1139
1140//=============================================================================================================
1141
1142MneSurfaceOld* MneSurfaceOrVolume::mne_read_bem_surface2(char *name, int which, int add_geometry, float *sigmap)
1143{
1144 return read_bem_surface(name,which,add_geometry,sigmap,FALSE);
1145}
1146
1147//=============================================================================================================
1148
1149MneSurfaceOld* MneSurfaceOrVolume::read_bem_surface(const QString &name, int which, int add_geometry, float *sigmap, bool check_too_many_neighbors)
1150/*
1151 * Read a Neuromag-style BEM surface description
1152 */
1153{
1154 QFile file(name);
1155 FiffStream::SPtr stream(new FiffStream(&file));
1156
1157 QList<FiffDirNode::SPtr> surfs;
1158 QList<FiffDirNode::SPtr> bems;
1159 FiffDirNode::SPtr node;
1160 FiffTag::SPtr t_pTag;
1161
1162 int id = -1;
1163 float **nodes = NULL;
1164 float **node_normals = NULL;
1165 int **triangles = NULL;
1166 int nnode,ntri;
1167 MneSurfaceOld* s = NULL;
1168 int k;
1169 int coord_frame = FIFFV_COORD_MRI;
1170 float sigma = -1.0;
1171 MatrixXf tmp_nodes;
1172 MatrixXi tmp_triangles;
1173
1174 if(!stream->open())
1175 goto bad;
1176 /*
1177 * Check for the existence of BEM coord frame
1178 */
1179 bems = stream->dirtree()->dir_tree_find(FIFFB_BEM);
1180 if (bems.size() > 0) {
1181 node = bems[0];
1182 if (node->find_tag(stream, FIFF_BEM_COORD_FRAME, t_pTag)) {
1183 coord_frame = *t_pTag->toInt();
1184 }
1185 }
1186 surfs = stream->dirtree()->dir_tree_find(FIFFB_BEM_SURF);
1187 if (surfs.size() == 0) {
1188 printf ("No BEM surfaces found in %s",name.toUtf8().constData());
1189 goto bad;
1190 }
1191 if (which >= 0) {
1192 for (k = 0; k < surfs.size(); ++k) {
1193 node = surfs[k];
1194 /*
1195 * Read the data from this node
1196 */
1197 if (node->find_tag(stream, FIFF_BEM_SURF_ID, t_pTag)) {
1198 id = *t_pTag->toInt();
1199 if (id == which)
1200 break;
1201 }
1202 }
1203 if (id != which) {
1204 printf("Desired surface not found in %s",name.toUtf8().constData());
1205 goto bad;
1206 }
1207 }
1208 else
1209 node = surfs[0];
1210 /*
1211 * Get the compulsory tags
1212 */
1213 if (!node->find_tag(stream, FIFF_BEM_SURF_NNODE, t_pTag))
1214 goto bad;
1215 nnode = *t_pTag->toInt();
1216
1217 if (!node->find_tag(stream, FIFF_BEM_SURF_NTRI, t_pTag))
1218 goto bad;
1219 ntri = *t_pTag->toInt();
1220
1221 if (!node->find_tag(stream, FIFF_BEM_SURF_NODES, t_pTag))
1222 goto bad;
1223 tmp_nodes = t_pTag->toFloatMatrix().transpose();
1224 nodes = ALLOC_CMATRIX_17(tmp_nodes.rows(),tmp_nodes.cols());
1225 fromFloatEigenMatrix_17(tmp_nodes, nodes);
1226
1227 if (node->find_tag(stream, FIFF_BEM_SURF_NORMALS, t_pTag)) {\
1228 MatrixXf tmp_node_normals = t_pTag->toFloatMatrix().transpose();
1229 node_normals = ALLOC_CMATRIX_17(tmp_node_normals.rows(),tmp_node_normals.cols());
1230 fromFloatEigenMatrix_17(tmp_node_normals, node_normals);
1231 }
1232
1233 if (!node->find_tag(stream, FIFF_BEM_SURF_TRIANGLES, t_pTag))
1234 goto bad;
1235 tmp_triangles = t_pTag->toIntMatrix().transpose();
1236 triangles = (int **)malloc(tmp_triangles.rows() * sizeof(int *));
1237 for (int i = 0; i < tmp_triangles.rows(); ++i)
1238 triangles[i] = (int *)malloc(tmp_triangles.cols() * sizeof(int));
1239 fromIntEigenMatrix_17(tmp_triangles, triangles);
1240
1241 if (node->find_tag(stream, FIFF_MNE_COORD_FRAME, t_pTag)) {
1242 coord_frame = *t_pTag->toInt();
1243 }
1244 else if (node->find_tag(stream, FIFF_BEM_COORD_FRAME, t_pTag)) {
1245 coord_frame = *t_pTag->toInt();
1246 }
1247 if (node->find_tag(stream, FIFF_BEM_SIGMA, t_pTag)) {
1248 sigma = *t_pTag->toFloat();
1249 }
1250
1251 stream->close();
1252
1253 s = (MneSurfaceOld*)mne_new_source_space(0);
1254 for (k = 0; k < ntri; k++) {
1255 triangles[k][0]--;
1256 triangles[k][1]--;
1257 triangles[k][2]--;
1258 }
1259 s->itris = triangles;
1260 s->id = which;
1261 s->coord_frame = coord_frame;
1262 s->rr = nodes; nodes = NULL;
1263 s->nn = node_normals; node_normals = NULL;
1264 s->ntri = ntri;
1265 s->np = nnode;
1266 s->curv = NULL;
1267 s->val = NULL;
1268
1269 if (add_geometry) {
1270 if (check_too_many_neighbors) {
1271 if (mne_source_space_add_geometry_info((MneSourceSpaceOld*)s,!s->nn) != OK)
1272 goto bad;
1273 }
1274 else {
1275 if (mne_source_space_add_geometry_info2((MneSourceSpaceOld*)s,!s->nn) != OK)
1276 goto bad;
1277 }
1278 }
1279 else if (s->nn == NULL) { /* Normals only */
1280 if (mne_add_vertex_normals((MneSourceSpaceOld*)s) != OK)
1281 goto bad;
1282 }
1283 else
1284 mne_add_triangle_data((MneSourceSpaceOld*)s);
1285
1286 s->nuse = s->np;
1287 s->inuse = MALLOC_17(s->np,int);
1288 s->vertno = MALLOC_17(s->np,int);
1289 for (k = 0; k < s->np; k++) {
1290 s->inuse[k] = TRUE;
1291 s->vertno[k] = k;
1292 }
1293 if (sigmap)
1294 *sigmap = sigma;
1295
1296 return s;
1297
1298bad : {
1299 FREE_CMATRIX_17(nodes);
1300 FREE_CMATRIX_17(node_normals);
1301 FREE_ICMATRIX_17(triangles);
1302 stream->close();
1303 return NULL;
1304 }
1305}
1306
1307//=============================================================================================================
1308
1309void MneSurfaceOrVolume::mne_triangle_coords(float *r, MneSurfaceOld* s, int tri, float *x, float *y, float *z)
1310/*
1311 * Compute the coordinates of a point within a triangle
1312 */
1313{
1314 double rr[3]; /* Vector from triangle corner #1 to r */
1315 double a,b,c,v1,v2,det;
1316 MneTriangle* this_tri;
1317
1318 this_tri = s->tris+tri;
1319
1320 VEC_DIFF_17(this_tri->r1,r,rr);
1321 *z = VEC_DOT_17(rr,this_tri->nn);
1322
1323 a = VEC_DOT_17(this_tri->r12,this_tri->r12);
1324 b = VEC_DOT_17(this_tri->r13,this_tri->r13);
1325 c = VEC_DOT_17(this_tri->r12,this_tri->r13);
1326
1327 v1 = VEC_DOT_17(rr,this_tri->r12);
1328 v2 = VEC_DOT_17(rr,this_tri->r13);
1329
1330 det = a*b - c*c;
1331
1332 *x = (b*v1 - c*v2)/det;
1333 *y = (a*v2 - c*v1)/det;
1334
1335 return;
1336}
1337
1338//=============================================================================================================
1339
1340int MneSurfaceOrVolume::nearest_triangle_point(float *r, MneSurfaceOld* s, void *user, int tri, float *x, float *y, float *z)
1341/*
1342 * Find the nearest point from a triangle
1343 */
1344{
1345
1346 double p,q,p0,q0,t0;
1347 double rr[3]; /* Vector from triangle corner #1 to r */
1348 double a,b,c,v1,v2,det;
1349 double best,dist,dist0;
1350 MneProjData* pd = (MneProjData*)user;
1351 MneTriangle* this_tri;
1352
1353 this_tri = s->tris+tri;
1354 VEC_DIFF_17(this_tri->r1,r,rr);
1355 dist = VEC_DOT_17(rr,this_tri->nn);
1356
1357 if (pd) {
1358 if (!pd->act[tri])
1359 return FALSE;
1360 a = pd->a[tri];
1361 b = pd->b[tri];
1362 c = pd->c[tri];
1363 }
1364 else {
1365 a = VEC_DOT_17(this_tri->r12,this_tri->r12);
1366 b = VEC_DOT_17(this_tri->r13,this_tri->r13);
1367 c = VEC_DOT_17(this_tri->r12,this_tri->r13);
1368 }
1369
1370 v1 = VEC_DOT_17(rr,this_tri->r12);
1371 v2 = VEC_DOT_17(rr,this_tri->r13);
1372
1373 det = a*b - c*c;
1374
1375 p = (b*v1 - c*v2)/det;
1376 q = (a*v2 - c*v1)/det;
1377 /*
1378 * If the point projects into the triangle we are done
1379 */
1380 if (p >= 0.0 && p <= 1.0 &&
1381 q >= 0.0 && q <= 1.0 &&
1382 q <= 1.0 - p) {
1383 *x = p;
1384 *y = q;
1385 *z = dist;
1386 return TRUE;
1387 }
1388 /*
1389 * Tough: must investigate the sides
1390 * We might do something intelligent here. However, for now it is ok
1391 * to do it in the hard way
1392 */
1393 /*
1394 * Side 1 -> 2
1395 */
1396 p0 = p + 0.5*(q * c)/a;
1397 if (p0 < 0.0)
1398 p0 = 0.0;
1399 else if (p0 > 1.0)
1400 p0 = 1.0;
1401 q0 = 0.0;
1402 dist0 = sqrt((p-p0)*(p-p0)*a +
1403 (q-q0)*(q-q0)*b +
1404 (p-p0)*(q-q0)*c +
1405 dist*dist);
1406 best = dist0;
1407 *x = p0;
1408 *y = q0;
1409 *z = dist0;
1410 /*
1411 * Side 2 -> 3
1412 */
1413 t0 = 0.5*((2.0*a-c)*(1.0-p) + (2.0*b-c)*q)/(a+b-c);
1414 if (t0 < 0.0)
1415 t0 = 0.0;
1416 else if (t0 > 1.0)
1417 t0 = 1.0;
1418 p0 = 1.0 - t0;
1419 q0 = t0;
1420 dist0 = sqrt((p-p0)*(p-p0)*a +
1421 (q-q0)*(q-q0)*b +
1422 (p-p0)*(q-q0)*c +
1423 dist*dist);
1424 if (dist0 < best) {
1425 best = dist0;
1426 *x = p0;
1427 *y = q0;
1428 *z = dist0;
1429 }
1430 /*
1431 * Side 1 -> 3
1432 */
1433 p0 = 0.0;
1434 q0 = q + 0.5*(p * c)/b;
1435 if (q0 < 0.0)
1436 q0 = 0.0;
1437 else if (q0 > 1.0)
1438 q0 = 1.0;
1439 dist0 = sqrt((p-p0)*(p-p0)*a +
1440 (q-q0)*(q-q0)*b +
1441 (p-p0)*(q-q0)*c +
1442 dist*dist);
1443 if (dist0 < best) {
1444 best = dist0;
1445 *x = p0;
1446 *y = q0;
1447 *z = dist0;
1448 }
1449 return TRUE;
1450}
1451
1452//=============================================================================================================
1453
1454void MneSurfaceOrVolume::project_to_triangle(MneSurfaceOld* s, int tri, float p, float q, float *r)
1455{
1456 int k;
1457 MneTriangle* this_tri;
1458
1459 this_tri = s->tris+tri;
1460
1461 for (k = 0; k < 3; k++)
1462 r[k] = this_tri->r1[k] + p*this_tri->r12[k] + q*this_tri->r13[k];
1463
1464 return;
1465}
1466
1467//=============================================================================================================
1468
1469int MneSurfaceOrVolume::mne_nearest_triangle_point(float *r, MneSurfaceOld* s, int tri, float *x, float *y, float *z)
1470/*
1471 * This is for external use
1472 */
1473{
1474 return nearest_triangle_point(r,s,NULL,tri,x,y,z);
1475}
1476
1477//=============================================================================================================
1478
1479int MneSurfaceOrVolume::mne_project_to_surface(MneSurfaceOld* s, void *proj_data, float *r, int project_it, float *distp)
1480/*
1481 * Project the point onto the closest point on the surface
1482 */
1483{
1484 float dist; /* Distance to the triangle */
1485 float p,q; /* Coordinates on the triangle */
1486 float p0,q0,dist0;
1487 int best;
1488 int k;
1489
1490 p0 = q0 = 0.0;
1491 dist0 = 0.0;
1492 for (best = -1, k = 0; k < s->ntri; k++) {
1493 if (nearest_triangle_point(r,s,proj_data,k,&p,&q,&dist)) {
1494 if (best < 0 || std::fabs(dist) < std::fabs(dist0)) {
1495 dist0 = dist;
1496 best = k;
1497 p0 = p;
1498 q0 = q;
1499 }
1500 }
1501 }
1502 if (best >= 0 && project_it)
1503 project_to_triangle(s,best,p0,q0,r);
1504 if (distp)
1505 *distp = dist0;
1506 return best;
1507}
1508
1509//=============================================================================================================
1510
1511void MneSurfaceOrVolume::mne_project_to_triangle(MneSurfaceOld* s,
1512 int best,
1513 float *r,
1514 float *proj)
1515/*
1516 * Project to a triangle provided that we know the best match already
1517 */
1518{
1519 float p,q,dist;
1520
1521 mne_nearest_triangle_point(r,s,best,&p,&q,&dist);
1522 project_to_triangle(s,best,p,q,proj);
1523
1524 return;
1525}
1526
1527//=============================================================================================================
1528
1529void MneSurfaceOrVolume::mne_find_closest_on_surface_approx(MneSurfaceOld* s, float **r, int np, int *nearest, float *dist, int nstep)
1530/*
1531 * Find the closest triangle on the surface for each point and the distance to it
1532 * This uses the values in nearest as approximations of the closest triangle
1533 */
1534{
1535 MneProjData* p = new MneProjData(s);
1536 int k,was;
1537 float mydist;
1538
1539 printf("%s for %d points %d steps...",nearest[0] < 0 ? "Closest" : "Approx closest",np,nstep);
1540
1541 for (k = 0; k < np; k++) {
1542 was = nearest[k];
1543 decide_search_restriction(s,p,nearest[k],nstep,r[k]);
1544 nearest[k] = mne_project_to_surface(s,p,r[k],0,dist ? dist+k : &mydist);
1545 if (nearest[k] < 0) {
1546 decide_search_restriction(s,p,-1,nstep,r[k]);
1547 nearest[k] = mne_project_to_surface(s,p,r[k],0,dist ? dist+k : &mydist);
1548 }
1549 }
1550 (void)was; // squash compiler warning, set but unused
1551
1552 printf("[done]\n");
1553 delete p;
1554 return;
1555}
1556
1557//=============================================================================================================
1558
1559void MneSurfaceOrVolume::decide_search_restriction(MneSurfaceOld* s,
1560 MneProjData* p,
1561 int approx_best, /* We know the best triangle approximately
1562 * already */
1563 int nstep,
1564 float *r)
1565/*
1566 * Restrict the search only to feasible triangles
1567 */
1568{
1569 int k;
1570 float diff[3],dist,mindist;
1571 int minvert;
1572
1573 for (k = 0; k < s->ntri; k++)
1574 p->act[k] = FALSE;
1575
1576 if (approx_best < 0) {
1577 /*
1578 * Search for the closest vertex
1579 */
1580 mindist = 1000.0;
1581 minvert = 0;
1582 for (k = 0; k < s->np; k++) {
1583 VEC_DIFF_17(r,s->rr[k],diff);
1584 dist = VEC_LEN_17(diff);
1585 if (dist < mindist && s->nneighbor_tri[k] > 0) {
1586 mindist = dist;
1587 minvert = k;
1588 }
1589 }
1590 }
1591 else {
1592 /*
1593 * Just use this triangle
1594 */
1595 MneTriangle* this_tri = NULL;
1596
1597 this_tri = s->tris+approx_best;
1598 VEC_DIFF_17(r,this_tri->r1,diff);
1599 mindist = VEC_LEN_17(diff);
1600 minvert = this_tri->vert[0];
1601
1602 VEC_DIFF_17(r,this_tri->r2,diff);
1603 dist = VEC_LEN_17(diff);
1604 if (dist < mindist) {
1605 mindist = dist;
1606 minvert = this_tri->vert[1];
1607 }
1608 VEC_DIFF_17(r,this_tri->r3,diff);
1609 dist = VEC_LEN_17(diff);
1610 if (dist < mindist) {
1611 mindist = dist;
1612 minvert = this_tri->vert[2];
1613 }
1614 }
1615 /*
1616 * Activate triangles in the neighborhood
1617 */
1618 activate_neighbors(s,minvert,p->act,nstep);
1619
1620 for (k = 0, p->nactive = 0; k < s->ntri; k++)
1621 if (p->act[k])
1622 p->nactive++;
1623 return;
1624}
1625
1626//=============================================================================================================
1627
1628void MneSurfaceOrVolume::activate_neighbors(MneSurfaceOld* s, int start, int *act, int nstep)
1629/*
1630 * Blessed recursion...
1631 */
1632{
1633 int k;
1634
1635 if (nstep == 0)
1636 return;
1637
1638 for (k = 0; k < s->nneighbor_tri[start]; k++)
1639 act[s->neighbor_tri[start][k]] = TRUE;
1640 for (k = 0; k < s->nneighbor_vert[start]; k++)
1641 activate_neighbors(s,s->neighbor_vert[start][k],act,nstep-1);
1642
1643 return;
1644}
1645
1646//=============================================================================================================
1647
1648int MneSurfaceOrVolume::mne_read_source_spaces(const QString &name, MneSourceSpaceOld* **spacesp, int *nspacep)
1649/*
1650 * Read source spaces from a FIFF file
1651 */
1652{
1653 QFile file(name);
1654 FiffStream::SPtr stream(new FiffStream(&file));
1655
1656 int nspace = 0;
1657 MneSourceSpaceOld* *spaces = NULL;
1658 MneSourceSpaceOld* new_space = NULL;
1659 QList<FiffDirNode::SPtr> sources;
1660 FiffDirNode::SPtr node;
1661 FiffTag::SPtr t_pTag;
1662 int j,k,p,q;
1663 int ntri;
1664 int *nearest = NULL;
1665 float *nearest_dist = NULL;
1666 int *nneighbors = NULL;
1667 int *neighbors = NULL;
1668 int *vol_dims = NULL;
1669
1670 if(!stream->open())
1671 goto bad;
1672
1673 sources = stream->dirtree()->dir_tree_find(FIFFB_MNE_SOURCE_SPACE);
1674 if (sources.size() == 0) {
1675 printf("No source spaces available here");
1676 goto bad;
1677 }
1678 for (j = 0; j < sources.size(); j++) {
1679 new_space = MneSurfaceOrVolume::mne_new_source_space(0);
1680 node = sources[j];
1681 /*
1682 * Get the mandatory data first
1683 */
1684 if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NPOINTS, t_pTag)) {
1685 goto bad;
1686 }
1687 new_space->np = *t_pTag->toInt();
1688 if (new_space->np == 0) {
1689 printf("No points in this source space");
1690 goto bad;
1691 }
1692 if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_POINTS, t_pTag)) {
1693 goto bad;
1694 }
1695 MatrixXf tmp_rr = t_pTag->toFloatMatrix().transpose();
1696 new_space->rr = ALLOC_CMATRIX_17(tmp_rr.rows(),tmp_rr.cols());
1697 fromFloatEigenMatrix_17(tmp_rr,new_space->rr);
1698 if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NORMALS, t_pTag)) {
1699 goto bad;
1700 }
1701 MatrixXf tmp_nn = t_pTag->toFloatMatrix().transpose();
1702 new_space->nn = ALLOC_CMATRIX_17(tmp_nn.rows(),tmp_nn.cols());
1703 fromFloatEigenMatrix_17(tmp_nn,new_space->nn);
1704 if (!node->find_tag(stream, FIFF_MNE_COORD_FRAME, t_pTag)) {
1705 new_space->coord_frame = FIFFV_COORD_MRI;
1706 }
1707 else {
1708 new_space->coord_frame = *t_pTag->toInt();
1709 }
1710 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_ID, t_pTag)) {
1711 new_space->id = *t_pTag->toInt();
1712 }
1713 if (node->find_tag(stream, FIFF_SUBJ_HIS_ID, t_pTag)) {
1714 new_space->subject = (char *)t_pTag->data();
1715 }
1716 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_TYPE, t_pTag)) {
1717 new_space->type = *t_pTag->toInt();
1718 }
1719 ntri = 0;
1720 if (node->find_tag(stream, FIFF_BEM_SURF_NTRI, t_pTag)) {
1721 ntri = *t_pTag->toInt();
1722 }
1723 else if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NTRI, t_pTag)) {
1724 ntri = *t_pTag->toInt();
1725 }
1726 if (ntri > 0) {
1727 int **itris = NULL;
1728
1729 if (!node->find_tag(stream, FIFF_BEM_SURF_TRIANGLES, t_pTag)) {
1730 if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_TRIANGLES, t_pTag)) {
1731 goto bad;
1732 }
1733 }
1734
1735 MatrixXi tmp_itris = t_pTag->toIntMatrix().transpose();
1736 itris = (int **)malloc(tmp_itris.rows() * sizeof(int *));
1737 for (int i = 0; i < tmp_itris.rows(); ++i)
1738 itris[i] = (int *)malloc(tmp_itris.cols() * sizeof(int));
1739 fromIntEigenMatrix_17(tmp_itris, itris);
1740
1741 for (p = 0; p < ntri; p++) { /* Adjust the numbering */
1742 itris[p][X_17]--;
1743 itris[p][Y_17]--;
1744 itris[p][Z_17]--;
1745 }
1746 new_space->itris = itris; itris = NULL;
1747 new_space->ntri = ntri;
1748 }
1749 if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NUSE, t_pTag)) {
1750 if (new_space->type == FIFFV_MNE_SPACE_VOLUME) {
1751 /*
1752 * Use all
1753 */
1754 new_space->nuse = new_space->np;
1755 new_space->inuse = MALLOC_17(new_space->nuse,int);
1756 new_space->vertno = MALLOC_17(new_space->nuse,int);
1757 for (k = 0; k < new_space->nuse; k++) {
1758 new_space->inuse[k] = TRUE;
1759 new_space->vertno[k] = k;
1760 }
1761 }
1762 else {
1763 /*
1764 * None in use
1765 * NOTE: The consequences of this change have to be evaluated carefully
1766 */
1767 new_space->nuse = 0;
1768 new_space->inuse = MALLOC_17(new_space->np,int);
1769 new_space->vertno = NULL;
1770 for (k = 0; k < new_space->np; k++)
1771 new_space->inuse[k] = FALSE;
1772 }
1773 }
1774 else {
1775 new_space->nuse = *t_pTag->toInt();
1776 if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_SELECTION, t_pTag)) {
1777 goto bad;
1778 }
1779
1780 qDebug() << "ToDo: Check whether new_space->inuse contains the right stuff!!! - use VectorXi instead";
1781 // new_space->inuse = t_pTag->toInt();
1782 new_space->inuse = MALLOC_17(new_space->np,int); //DEBUG
1783 if (new_space->nuse > 0) {
1784 new_space->vertno = MALLOC_17(new_space->nuse,int);
1785 for (k = 0, p = 0; k < new_space->np; k++) {
1786 new_space->inuse[k] = t_pTag->toInt()[k]; //DEBUG
1787 if (new_space->inuse[k])
1788 new_space->vertno[p++] = k;
1789 }
1790 }
1791 else {
1792 FREE_17(new_space->vertno);
1793 new_space->vertno = NULL;
1794 }
1795 /*
1796 * Selection triangulation
1797 */
1798 ntri = 0;
1799 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NUSE_TRI, t_pTag)) {
1800 ntri = *t_pTag->toInt();
1801 }
1802 if (ntri > 0) {
1803 int **itris = NULL;
1804
1805 if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_USE_TRIANGLES, t_pTag)) {
1806 goto bad;
1807 }
1808
1809 MatrixXi tmp_itris = t_pTag->toIntMatrix().transpose();
1810 itris = (int **)malloc(tmp_itris.rows() * sizeof(int *));
1811 for (int i = 0; i < tmp_itris.rows(); ++i)
1812 itris[i] = (int *)malloc(tmp_itris.cols() * sizeof(int));
1813 fromIntEigenMatrix_17(tmp_itris, itris);
1814 for (p = 0; p < ntri; p++) { /* Adjust the numbering */
1815 itris[p][X_17]--;
1816 itris[p][Y_17]--;
1817 itris[p][Z_17]--;
1818 }
1819 new_space->use_itris = itris; itris = NULL;
1820 new_space->nuse_tri = ntri;
1821 }
1822 /*
1823 * The patch information becomes relevant here
1824 */
1825 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NEAREST, t_pTag)) {
1826 nearest = t_pTag->toInt();
1827 new_space->nearest = MALLOC_17(new_space->np,MneNearest);
1828 for (k = 0; k < new_space->np; k++) {
1829 new_space->nearest[k].vert = k;
1830 new_space->nearest[k].nearest = nearest[k];
1831 new_space->nearest[k].patch = NULL;
1832 }
1833
1834 if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NEAREST_DIST, t_pTag)) {
1835 goto bad;
1836 }
1837 qDebug() << "ToDo: Check whether nearest_dist contains the right stuff!!! - use VectorXf instead";
1838 nearest_dist = t_pTag->toFloat();
1839 for (k = 0; k < new_space->np; k++) {
1840 new_space->nearest[k].dist = nearest_dist[k];
1841 }
1842 // FREE_17(nearest); nearest = NULL;
1843 // FREE_17(nearest_dist); nearest_dist = NULL;
1844 }
1845 /*
1846 * We may have the distance matrix
1847 */
1848 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_DIST_LIMIT, t_pTag)) {
1849 new_space->dist_limit = *t_pTag->toFloat();
1850 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_DIST, t_pTag)) {
1851 // SparseMatrix<double> tmpSparse = t_pTag->toSparseFloatMatrix();
1852 FiffSparseMatrix* dist = FiffSparseMatrix::fiff_get_float_sparse_matrix(t_pTag);
1853 new_space->dist = dist->mne_add_upper_triangle_rcs();
1854 delete dist;
1855 if (!new_space->dist) {
1856 goto bad;
1857 }
1858 }
1859 else
1860 new_space->dist_limit = 0.0;
1861 }
1862 }
1863 /*
1864 * For volume source spaces we might have the neighborhood information
1865 */
1866 if (new_space->type == FIFFV_MNE_SPACE_VOLUME) {
1867 int ntot,nvert,ntot_count,nneigh;
1868 int *neigh;
1869
1870 nneighbors = neighbors = NULL;
1871 ntot = nvert = 0;
1872 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NEIGHBORS, t_pTag)) {
1873 qDebug() << "ToDo: Check whether neighbors contains the right stuff!!! - use VectorXi instead";
1874 neighbors = t_pTag->toInt();
1875 ntot = t_pTag->size()/sizeof(fiff_int_t);
1876 }
1877 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NNEIGHBORS, t_pTag)) {
1878 qDebug() << "ToDo: Check whether nneighbors contains the right stuff!!! - use VectorXi instead";
1879 nneighbors = t_pTag->toInt();
1880 nvert = t_pTag->size()/sizeof(fiff_int_t);
1881 }
1882 if (neighbors && nneighbors) {
1883 if (nvert != new_space->np) {
1884 printf("Inconsistent neighborhood data in file.");
1885 goto bad;
1886 }
1887 for (k = 0, ntot_count = 0; k < nvert; k++)
1888 ntot_count += nneighbors[k];
1889 if (ntot_count != ntot) {
1890 printf("Inconsistent neighborhood data in file.");
1891 goto bad;
1892 }
1893 new_space->nneighbor_vert = MALLOC_17(nvert,int);
1894 new_space->neighbor_vert = MALLOC_17(nvert,int *);
1895 for (k = 0, q = 0; k < nvert; k++) {
1896 new_space->nneighbor_vert[k] = nneigh = nneighbors[k];
1897 new_space->neighbor_vert[k] = neigh = MALLOC_17(nneigh,int);
1898 for (p = 0; p < nneigh; p++,q++)
1899 neigh[p] = neighbors[q];
1900 }
1901 }
1902 FREE_17(neighbors);
1903 FREE_17(nneighbors);
1904 nneighbors = neighbors = NULL;
1905 /*
1906 * There might be a coordinate transformation and dimensions
1907 */
1908 new_space->voxel_surf_RAS_t = FiffCoordTransOld::mne_read_transform_from_node(stream, node, FIFFV_MNE_COORD_MRI_VOXEL, FIFFV_MNE_COORD_SURFACE_RAS);
1909 if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_VOXEL_DIMS, t_pTag)) {
1910 qDebug() << "ToDo: Check whether vol_dims contains the right stuff!!! - use VectorXi instead";
1911 vol_dims = t_pTag->toInt();
1912 }
1913 if (vol_dims)
1914 VEC_COPY_17(new_space->vol_dims,vol_dims);
1915 {
1916 QList<FiffDirNode::SPtr> mris = node->dir_tree_find(FIFFB_MNE_PARENT_MRI_FILE);
1917
1918 if (mris.size() == 0) { /* The old way */
1919 new_space->MRI_surf_RAS_RAS_t = FiffCoordTransOld::mne_read_transform_from_node(stream, node, FIFFV_MNE_COORD_SURFACE_RAS, FIFFV_MNE_COORD_RAS);
1920 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_MRI_FILE, t_pTag)) {
1921 qDebug() << "ToDo: Check whether new_space->MRI_volume contains the right stuff!!! - use QString instead";
1922 new_space->MRI_volume = (char *)t_pTag->data();
1923 }
1924 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_INTERPOLATOR, t_pTag)) {
1925 new_space->interpolator = FiffSparseMatrix::fiff_get_float_sparse_matrix(t_pTag);
1926 }
1927 }
1928 else {
1929 if (node->find_tag(stream, FIFF_MNE_FILE_NAME, t_pTag)) {
1930 new_space->MRI_volume = (char *)t_pTag->data();
1931 }
1932 new_space->MRI_surf_RAS_RAS_t = FiffCoordTransOld::mne_read_transform_from_node(stream, mris[0], FIFFV_MNE_COORD_SURFACE_RAS, FIFFV_MNE_COORD_RAS);
1933 new_space->MRI_voxel_surf_RAS_t = FiffCoordTransOld::mne_read_transform_from_node(stream, mris[0], FIFFV_MNE_COORD_MRI_VOXEL, FIFFV_MNE_COORD_SURFACE_RAS);
1934
1935 if (mris[0]->find_tag(stream, FIFF_MNE_SOURCE_SPACE_INTERPOLATOR, t_pTag)) {
1936 new_space->interpolator = FiffSparseMatrix::fiff_get_float_sparse_matrix(t_pTag);
1937 }
1938 if (mris[0]->find_tag(stream, FIFF_MRI_WIDTH, t_pTag)) {
1939 new_space->MRI_vol_dims[0] = *t_pTag->toInt();
1940 }
1941 if (mris[0]->find_tag(stream, FIFF_MRI_HEIGHT, t_pTag)) {
1942 new_space->MRI_vol_dims[1] = *t_pTag->toInt();
1943 }
1944 if (mris[0]->find_tag(stream, FIFF_MRI_DEPTH, t_pTag)) {
1945 new_space->MRI_vol_dims[2] = *t_pTag->toInt();
1946 }
1947 }
1948 }
1949 }
1950 mne_add_triangle_data(new_space);
1951 spaces = REALLOC_17(spaces,nspace+1,MneSourceSpaceOld*);
1952 spaces[nspace++] = new_space;
1953 new_space = NULL;
1954 }
1955 stream->close();
1956
1957 *spacesp = spaces;
1958 *nspacep = nspace;
1959
1960 return FIFF_OK;
1961
1962bad : {
1963 stream->close();
1964 delete new_space;
1965 for (k = 0; k < nspace; k++)
1966 delete spaces[k];
1967 FREE_17(spaces);
1968 FREE_17(nearest);
1969 FREE_17(nearest_dist);
1970 FREE_17(neighbors);
1971 FREE_17(nneighbors);
1972 FREE_17(vol_dims);
1973
1974 return FIFF_FAIL;
1975 }
1976}
1977
1978//=============================================================================================================
1979
1980void MneSurfaceOrVolume::mne_source_space_update_inuse(MneSourceSpaceOld* s, int *new_inuse)
1981/*
1982 * Update the active vertices
1983 */
1984{
1985 int k,p,nuse;
1986
1987 if (!s)
1988 return;
1989
1990 FREE_17(s->inuse); s->inuse = new_inuse;
1991
1992 for (k = 0, nuse = 0; k < s->np; k++)
1993 if (s->inuse[k])
1994 nuse++;
1995
1996 s->nuse = nuse;
1997 if (s->nuse > 0) {
1998 s->vertno = REALLOC_17(s->vertno,s->nuse,int);
1999 for (k = 0, p = 0; k < s->np; k++)
2000 if (s->inuse[k])
2001 s->vertno[p++] = k;
2002 }
2003 else {
2004 FREE_17(s->vertno);
2005 s->vertno = NULL;
2006 }
2007 return;
2008}
2009
2010//=============================================================================================================
2011
2012int MneSurfaceOrVolume::mne_is_left_hemi_source_space(MneSourceSpaceOld* s)
2013/*
2014 * Left or right hemisphere?
2015 */
2016{
2017 int k;
2018 float xave;
2019
2020 for (k = 0, xave = 0.0; k < s->np; k++)
2021 xave += s->rr[k][0];
2022 if (xave < 0.0)
2023 return TRUE;
2024 else
2025 return FALSE;
2026}
2027
2028//=============================================================================================================
2029
2030int MneSurfaceOrVolume::mne_transform_source_space(MneSourceSpaceOld* ss, FiffCoordTransOld *t)
2031/*
2032 * Transform source space data into another coordinate frame
2033 */
2034{
2035 int k;
2036 if (ss == NULL)
2037 return OK;
2038 if (ss->coord_frame == t->to)
2039 return OK;
2040 if (ss->coord_frame != t->from) {
2041 printf("Coordinate transformation does not match with the source space coordinate system.");
2042 return FAIL;
2043 }
2044 for (k = 0; k < ss->np; k++) {
2045 FiffCoordTransOld::fiff_coord_trans(ss->rr[k],t,FIFFV_MOVE);
2046 FiffCoordTransOld::fiff_coord_trans(ss->nn[k],t,FIFFV_NO_MOVE);
2047 }
2048 if (ss->tris) {
2049 for (k = 0; k < ss->ntri; k++)
2050 FiffCoordTransOld::fiff_coord_trans(ss->tris[k].nn,t,FIFFV_NO_MOVE);
2051 }
2052 ss->coord_frame = t->to;
2053 return OK;
2054}
2055
2056//=============================================================================================================
2057
2058int MneSurfaceOrVolume::mne_transform_source_spaces_to(int coord_frame, FiffCoordTransOld *t, MneSourceSpaceOld* *spaces, int nspace)
2059/*
2060 * Facilitate the transformation of the source spaces
2061 */
2062{
2063 MneSourceSpaceOld* s;
2064 int k;
2065 FiffCoordTransOld* my_t;
2066
2067 for (k = 0; k < nspace; k++) {
2068 s = spaces[k];
2069 if (s->coord_frame != coord_frame) {
2070 if (t) {
2071 if (s->coord_frame == t->from && t->to == coord_frame) {
2072 if (mne_transform_source_space(s,t) != OK)
2073 return FAIL;
2074 }
2075 else if (s->coord_frame == t->to && t->from == coord_frame) {
2076 my_t = t->fiff_invert_transform();
2077 if (mne_transform_source_space(s,my_t) != OK) {
2078 FREE_17(my_t);
2079 return FAIL;
2080 }
2081 FREE_17(my_t);
2082 }
2083 else {
2084 printf("Could not transform a source space because of transformation incompatibility.");
2085 return FAIL;
2086 }
2087 }
2088 else {
2089 printf("Could not transform a source space because of missing coordinate transformation.");
2090 return FAIL;
2091 }
2092 }
2093 }
2094 return OK;
2095}
2096
2097//=============================================================================================================
2098
2099void MneSurfaceOrVolume::enable_all_sources(MneSourceSpaceOld* s)
2100{
2101 int k;
2102 for (k = 0; k < s->np; k++)
2103 s->inuse[k] = TRUE;
2104 s->nuse = s->np;
2105 return;
2106}
2107
2108//=============================================================================================================
2109
2110#define LH_LABEL_TAG "-lh.label"
2111#define RH_LABEL_TAG "-rh.label"
2112
2113int MneSurfaceOrVolume::restrict_sources_to_labels(MneSourceSpaceOld* *spaces, int nspace, const QStringList& labels, int nlabel)
2114/*
2115 * Pick only sources within a label
2116 */
2117{
2118 MneSourceSpaceOld* lh = NULL;
2119 MneSourceSpaceOld* rh = NULL;
2120 MneSourceSpaceOld* sp;
2121 int *lh_inuse = NULL;
2122 int *rh_inuse = NULL;
2123 int *sel = NULL;
2124 int nsel;
2125 int *inuse;
2126 int k,p;
2127
2128 if (nlabel == 0)
2129 return OK;
2130
2131 for (k = 0; k < nspace; k++) {
2132 if (mne_is_left_hemi_source_space(spaces[k])) {
2133 lh = spaces[k];
2134 FREE_17(lh_inuse);
2135 lh_inuse = MALLOC_17(lh->np,int);
2136 for (p = 0; p < lh->np; p++)
2137 lh_inuse[p] = 0;
2138 }
2139 else {
2140 rh = spaces[k];
2141 FREE_17(rh_inuse);
2142 rh_inuse = MALLOC_17(rh->np,int);
2143 for (p = 0; p < rh->np; p++)
2144 rh_inuse[p] = 0;
2145 }
2146 }
2147 /*
2148 * Go through each label file
2149 */
2150 for (k = 0; k < nlabel; k++) {
2151 /*
2152 * Which hemi?
2153 */
2154 if (labels[k].contains(LH_LABEL_TAG)){ //strstr(labels[k],LH_LABEL_TAG) != NULL) {
2155 sp = lh;
2156 inuse = lh_inuse;
2157 }
2158 else if (labels[k].contains(RH_LABEL_TAG)){ //strstr(labels[k],RH_LABEL_TAG) != NULL) {
2159 sp = rh;
2160 inuse = rh_inuse;
2161 }
2162 else {
2163 printf("\tWarning: cannot assign label file %s to a hemisphere.\n",labels[k].toUtf8().constData());
2164 continue;
2165 }
2166 if (sp) {
2167 if (mne_read_label(labels[k],NULL,&sel,&nsel) == FAIL)
2168 goto bad;
2169 for (p = 0; p < nsel; p++) {
2170 if (sel[p] >= 0 && sel[p] < sp->np)
2171 inuse[sel[p]] = sp->inuse[sel[p]];
2172 else
2173 printf("vertex number out of range in %s (%d vs %d)\n",
2174 labels[k].toUtf8().constData(),sel[p],sp->np);
2175 }
2176 FREE_17(sel); sel = NULL;
2177 printf("Processed label file %s\n",labels[k].toUtf8().constData());
2178 }
2179 }
2180 mne_source_space_update_inuse(lh,lh_inuse);
2181 mne_source_space_update_inuse(rh,rh_inuse);
2182 return OK;
2183
2184bad : {
2185 FREE_17(lh_inuse);
2186 FREE_17(rh_inuse);
2187 FREE_17(sel);
2188 return FAIL;
2189 }
2190}
2191
2192//=============================================================================================================
2193
2194int MneSurfaceOrVolume::mne_find_sources_in_label(char *label, MneSourceSpaceOld* s, int off, int **selp, int *nselp) /* How many selected? */
2195/*
2196 * Find the source points within a label
2197 */
2198{
2199 FILE *in = NULL;
2200 int res = FAIL;
2201
2202 int nsel = 0;
2203 int *sel = NULL;
2204
2205 int k,p,pp,nlabel,q;
2206 char c;
2207 float fdum;
2208 /*
2209 * Read the label file
2210 */
2211 if ((in = fopen(label,"r")) == NULL) {
2212 qCritical("%s", label);//err_set_sys_error(label);
2213 goto out;
2214 }
2215 c = fgetc(in);
2216 if (c !='#') {
2217 qCritical("Label file does not start correctly.");
2218 goto out;
2219 }
2220 for (c = fgetc(in); c != '\n'; c = fgetc(in))
2221 ;
2222 if (fscanf(in,"%d",&nlabel) != 1) {
2223 qCritical("Could not read the number of labelled points.");
2224 goto out;
2225 }
2226#ifdef DEBUG
2227 printf("\t%d points in label %s\n",nlabel,label);
2228#endif
2229 for (k = 0; k < nlabel; k++) {
2230 if (fscanf(in,"%d %g %g %g %g",&p,
2231 &fdum, &fdum, &fdum, &fdum) != 5) {
2232 qCritical("Could not read label point # %d",k+1);
2233 goto out;
2234 }
2235 if (p < 0 || p >= s->np) {
2236 qCritical("Source index out of range %d (range 0..%d)\n",p,s->np-1);
2237 goto out;
2238 }
2239 if (s->inuse[p]) {
2240 for (pp = 0, q = 0; pp < p; pp++) {
2241 if (s->inuse[pp])
2242 q++;
2243 }
2244 sel = REALLOC_17(sel,nsel+1,int);
2245 sel[nsel++] = q + off;
2246 }
2247 }
2248 *nselp = nsel;
2249 *selp = sel;
2250 res = OK;
2251
2252out : {
2253 if (in)
2254 fclose(in);
2255 if (res != OK) {
2256 FREE_17(sel);
2257 *selp = NULL;
2258 *nselp = 0;
2259 }
2260 return res;
2261 }
2262}
2263
2264//=============================================================================================================
2265
2266int MneSurfaceOrVolume::mne_read_label(const QString& label, char **commentp, int **selp, int *nselp) /* How many? */
2267/*
2268 * Find the source points within a label
2269 */
2270{
2271 FILE *in = NULL;
2272 int res = FAIL;
2273
2274 int nsel = 0;
2275 int *sel = NULL;
2276
2277 int k,p,nlabel;
2278 char c;
2279 char *comment = NULL;
2280 float fdum;
2281 /*
2282 * Read the label file
2283 */
2284 if ((in = fopen(label.toUtf8().constData(),"r")) == NULL) {
2285 qCritical() << label;//err_set_sys_error(label);
2286 goto out;
2287 }
2288 for (k = 0; k < 2; k++) {
2289 rewind(in);
2290 c = fgetc(in);
2291 if (c !='#') {
2292 qCritical("Label file does not start correctly.");
2293 goto out;
2294 }
2295 for (c = fgetc(in); c == ' ' && c != '\n'; c = fgetc(in))
2296 ;
2297 ungetc(c,in);
2298 if (k == 0) {
2299 for (c = fgetc(in), p = 0; c != '\n'; c = fgetc(in), p++)
2300 ;
2301 }
2302 else {
2303 for (c = fgetc(in); c != '\n'; c = fgetc(in))
2304 *comment++ = c;
2305 *comment = '\0';
2306 }
2307 if (!commentp)
2308 break;
2309 if (p == 0) {
2310 *commentp = NULL;
2311 break;
2312 }
2313 if (k == 0)
2314 *commentp = comment = MALLOC_17(p+1,char);
2315 }
2316 if (fscanf(in,"%d",&nlabel) != 1) {
2317 qCritical("Could not read the number of labelled points.");
2318 goto out;
2319 }
2320 for (k = 0; k < nlabel; k++) {
2321 if (fscanf(in,"%d %g %g %g %g",&p,
2322 &fdum, &fdum, &fdum, &fdum) != 5) {
2323 qCritical("Could not read label point # %d",k+1);
2324 goto out;
2325 }
2326 sel = REALLOC_17(sel,nsel+1,int);
2327 sel[nsel++] = p;
2328 }
2329 *nselp = nsel;
2330 *selp = sel;
2331 res = OK;
2332
2333out : {
2334 if (in)
2335 fclose(in);
2336 if (res != OK) {
2337 FREE_17(sel);
2338 *selp = NULL;
2339 *nselp = 0;
2340 }
2341 return res;
2342 }
2343}
2344
2345//=============================================================================================================
2346
2347int MneSurfaceOrVolume::mne_write_label(char *label, char *comment, int *sel, int nsel, float **rr) /* Locations of the nodes in MRI coords */
2348/*
2349 * Find the source points within a label
2350 */
2351{
2352 FILE *out = NULL;
2353 int res = FAIL;
2354 int k;
2355 float fdum = 0.0;
2356 /*
2357 * Read the label file
2358 */
2359 if ((out = fopen(label,"w")) == NULL) {
2360 qCritical("%s", label);//err_set_sys_error(label);
2361 goto out;
2362 }
2363 if (comment == NULL)
2364 fprintf(out,"# Label file created by the MNE software\n");
2365 else
2366 fprintf(out,"# %s\n",comment);
2367
2368 fprintf(out,"%d\n",nsel);
2369 if (rr != NULL)
2370 for (k = 0; k < nsel; k++)
2371 fprintf(out,"%d %.2f %.2f %.2f %g\n",sel[k],
2372 1000*rr[sel[k]][0],
2373 1000*rr[sel[k]][1],
2374 1000*rr[sel[k]][2],fdum);
2375 else
2376 for (k = 0; k < nsel; k++)
2377 fprintf(out,"%d %.2f %.2f %.2f %g\n",sel[k],fdum,fdum,fdum,fdum);
2378 res = OK;
2379
2380out : {
2381 if (out)
2382 fclose(out);
2383 if (res != OK)
2384 unlink(label);
2385 return res;
2386 }
2387}
2388
2389//=============================================================================================================
2390
2391void MneSurfaceOrVolume::mne_add_triangle_data(MneSourceSpaceOld* s)
2392/*
2393 * Add the triangle data structures
2394 */
2395{
2396 int k;
2397 MneTriangle* tri;
2398
2399 if (!s || s->type != MNE_SOURCE_SPACE_SURFACE)
2400 return;
2401
2402 FREE_17(s->tris); s->tris = NULL;
2403 FREE_17(s->use_tris); s->use_tris = NULL;
2404 /*
2405 * Add information for the complete triangulation
2406 */
2407 if (s->itris && s->ntri > 0) {
2408 s->tris = MALLOC_17(s->ntri,MneTriangle);
2409 s->tot_area = 0.0;
2410 for (k = 0, tri = s->tris; k < s->ntri; k++, tri++) {
2411 tri->vert = s->itris[k];
2412 tri->r1 = s->rr[tri->vert[0]];
2413 tri->r2 = s->rr[tri->vert[1]];
2414 tri->r3 = s->rr[tri->vert[2]];
2415 MneTriangle::add_triangle_data(tri);
2416 s->tot_area += tri->area;
2417 }
2418#ifdef TRIANGLE_SIZE_WARNING
2419 for (k = 0, tri = s->tris; k < s->ntri; k++, tri++)
2420 if (tri->area < 1e-5*s->tot_area/s->ntri)
2421 printf("Warning: Triangle area is only %g um^2 (%.5f %% of expected average)\n",
2422 1e12*tri->area,100*s->ntri*tri->area/s->tot_area);
2423#endif
2424 }
2425#ifdef DEBUG
2426 printf("\ttotal area = %-.1f cm^2\n",1e4*s->tot_area);
2427#endif
2428 /*
2429 * Add information for the selected subset if applicable
2430 */
2431 if (s->use_itris && s->nuse_tri > 0) {
2432 s->use_tris = MALLOC_17(s->nuse_tri,MneTriangle);
2433 for (k = 0, tri = s->use_tris; k < s->nuse_tri; k++, tri++) {
2434 tri->vert = s->use_itris[k];
2435 tri->r1 = s->rr[tri->vert[0]];
2436 tri->r2 = s->rr[tri->vert[1]];
2437 tri->r3 = s->rr[tri->vert[2]];
2438 MneTriangle::add_triangle_data(tri);
2439 }
2440 }
2441 return;
2442}
2443
2444//=============================================================================================================
2445
2446void MneSurfaceOrVolume::mne_compute_cm(float **rr, int np, float *cm)
2447/*
2448 * Compute the center of mass of a set of points
2449 */
2450{
2451 int q;
2452 cm[0] = cm[1] = cm[2] = 0.0;
2453 for (q = 0; q < np; q++) {
2454 cm[0] += rr[q][0];
2455 cm[1] += rr[q][1];
2456 cm[2] += rr[q][2];
2457 }
2458 if (np > 0) {
2459 cm[0] = cm[0]/np;
2460 cm[1] = cm[1]/np;
2461 cm[2] = cm[2]/np;
2462 }
2463 return;
2464}
2465
2466//=============================================================================================================
2467
2468void MneSurfaceOrVolume::mne_compute_surface_cm(MneSurfaceOld *s)
2469/*
2470 * Compute the center of mass of a surface
2471 */
2472{
2473 if (!s)
2474 return;
2475
2476 mne_compute_cm(s->rr,s->np,s->cm);
2477 return;
2478}
2479
2480//=============================================================================================================
2481
2482void MneSurfaceOrVolume::calculate_vertex_distances(MneSourceSpaceOld* s)
2483{
2484 int k,p,ndist;
2485 float *dist,diff[3];
2486 int *neigh, nneigh;
2487
2488 if (!s->neighbor_vert || !s->nneighbor_vert)
2489 return;
2490
2491 if (s->vert_dist) {
2492 for (k = 0; k < s->np; k++)
2493 FREE_17(s->vert_dist[k]);
2494 FREE_17(s->vert_dist);
2495 }
2496 s->vert_dist = MALLOC_17(s->np,float *);
2497 printf("\tDistances between neighboring vertices...");
2498 for (k = 0, ndist = 0; k < s->np; k++) {
2499 s->vert_dist[k] = dist = MALLOC_17(s->nneighbor_vert[k],float);
2500 neigh = s->neighbor_vert[k];
2501 nneigh = s->nneighbor_vert[k];
2502 for (p = 0; p < nneigh; p++) {
2503 if (neigh[p] >= 0) {
2504 VEC_DIFF_17(s->rr[k],s->rr[neigh[p]],diff);
2505 dist[p] = VEC_LEN_17(diff);
2506 }
2507 else
2508 dist[p] = -1.0;
2509 ndist++;
2510 }
2511 }
2512 printf("[%d distances done]\n",ndist);
2513 return;
2514}
2515
2516//=============================================================================================================
2517
2518int MneSurfaceOrVolume::mne_add_vertex_normals(MneSourceSpaceOld* s)
2519{
2520 int k,c,p;
2521 int *ii;
2522 float w,size;
2523 MneTriangle* tri;
2524
2525 if (!s || s->type != MNE_SOURCE_SPACE_SURFACE)
2526 return OK;
2527 /*
2528 * Reallocate the stuff and initialize
2529 */
2530 FREE_CMATRIX_17(s->nn);
2531 s->nn = ALLOC_CMATRIX_17(s->np,3);
2532
2533 for (k = 0; k < s->np; k++) {
2534 s->nn[k][X_17] = s->nn[k][Y_17] = s->nn[k][Z_17] = 0.0;
2535 }
2536 /*
2537 * One pass through the triangles will do it
2538 */
2539 MneSurfaceOrVolume::mne_add_triangle_data(s);
2540 for (p = 0, tri = s->tris; p < s->ntri; p++, tri++) {
2541 ii = tri->vert;
2542 w = 1.0; /* This should be related to the triangle size */
2543 /*
2544 * Then the vertex normals
2545 */
2546 for (k = 0; k < 3; k++)
2547 for (c = 0; c < 3; c++)
2548 s->nn[ii[k]][c] += w*tri->nn[c];
2549 }
2550 for (k = 0; k < s->np; k++) {
2551 size = VEC_LEN_17(s->nn[k]);
2552 if (size > 0.0)
2553 for (c = 0; c < 3; c++)
2554 s->nn[k][c] = s->nn[k][c]/size;
2555 }
2556 mne_compute_surface_cm((MneSurfaceOld*)s);
2557 return OK;
2558}
2559
2560//=============================================================================================================
2561
2562int MneSurfaceOrVolume::add_geometry_info(MneSourceSpaceOld* s, int do_normals, int *border, int check_too_many_neighbors)
2563/*
2564 * Add vertex normals and neighbourhood information
2565 */
2566{
2567 int k,c,p,q;
2568 int vert;
2569 int *ii;
2570 int *neighbors,nneighbors;
2571 float w,size;
2572 int found;
2573 int nfix_distinct,nfix_no_neighbors,nfix_defect;
2574 MneTriangle* tri;
2575
2576 if (!s)
2577 return OK;
2578
2579 if (s->type == MNE_SOURCE_SPACE_VOLUME) {
2580 calculate_vertex_distances(s);
2581 return OK;
2582 }
2583 if (s->type != MNE_SOURCE_SPACE_SURFACE)
2584 return OK;
2585 /*
2586 * Reallocate the stuff and initialize
2587 */
2588 if (do_normals) {
2589 FREE_CMATRIX_17(s->nn);
2590 s->nn = ALLOC_CMATRIX_17(s->np,3);
2591 }
2592 if (s->neighbor_tri) {
2593 for (k = 0; k < s->np; k++)
2594 FREE_17(s->neighbor_tri[k]);
2595 FREE_17(s->neighbor_tri);
2596 }
2597 FREE_17(s->nneighbor_tri);
2598 s->neighbor_tri = MALLOC_17(s->np,int *);
2599 s->nneighbor_tri = MALLOC_17(s->np,int);
2600
2601 for (k = 0; k < s->np; k++) {
2602 s->neighbor_tri[k] = NULL;
2603 s->nneighbor_tri[k] = 0;
2604 if (do_normals)
2605 s->nn[k][X_17] = s->nn[k][Y_17] = s->nn[k][Z_17] = 0.0;
2606 }
2607 /*
2608 * One pass through the triangles will do it
2609 */
2610 mne_add_triangle_data(s);
2611 for (p = 0, tri = s->tris; p < s->ntri; p++, tri++)
2612 if (tri->area == 0)
2613 printf("\tWarning : zero size triangle # %d\n",p);
2614 printf("\tTriangle ");
2615 if (do_normals)
2616 printf("and vertex ");
2617 printf("normals and neighboring triangles...");
2618 for (p = 0, tri = s->tris; p < s->ntri; p++, tri++) {
2619 ii = tri->vert;
2620 w = 1.0; /* This should be related to the triangle size */
2621 for (k = 0; k < 3; k++) {
2622 /*
2623 * Then the vertex normals
2624 */
2625 if (do_normals)
2626 for (c = 0; c < 3; c++)
2627 s->nn[ii[k]][c] += w*tri->nn[c];
2628 /*
2629 * Add to the list of neighbors
2630 */
2631 s->neighbor_tri[ii[k]] = REALLOC_17(s->neighbor_tri[ii[k]],
2632 s->nneighbor_tri[ii[k]]+1,int);
2633 s->neighbor_tri[ii[k]][s->nneighbor_tri[ii[k]]] = p;
2634 s->nneighbor_tri[ii[k]]++;
2635 }
2636 }
2637 nfix_no_neighbors = 0;
2638 nfix_defect = 0;
2639 for (k = 0; k < s->np; k++) {
2640 if (s->nneighbor_tri[k] <= 0) {
2641 if (!border || !border[k]) {
2642#ifdef STRICT_ERROR
2643 err_printf_set_error("Vertex %d does not have any neighboring triangles!",k);
2644 return FAIL;
2645#else
2646#ifdef REPORT_WARNINGS
2647 printf("Warning: Vertex %d does not have any neighboring triangles!\n",k);
2648#endif
2649#endif
2650 nfix_no_neighbors++;
2651 }
2652 }
2653 else if (s->nneighbor_tri[k] < 3 && !border) {
2654#ifdef REPORT_WARNINGS
2655 printf("\n\tTopological defect: Vertex %d has only %d neighboring triangle%s Vertex omitted.\n\t",
2656 k,s->nneighbor_tri[k],s->nneighbor_tri[k] > 1 ? "s." : ".");
2657#endif
2658 nfix_defect++;
2659 s->nneighbor_tri[k] = 0;
2660 FREE_17(s->neighbor_tri[k]);
2661 s->neighbor_tri[k] = NULL;
2662 }
2663 }
2664 /*
2665 * Scale the vertex normals to unit length
2666 */
2667 for (k = 0; k < s->np; k++)
2668 if (s->nneighbor_tri[k] > 0) {
2669 size = VEC_LEN_17(s->nn[k]);
2670 if (size > 0.0)
2671 for (c = 0; c < 3; c++)
2672 s->nn[k][c] = s->nn[k][c]/size;
2673 }
2674 printf("[done]\n");
2675 /*
2676 * Determine the neighboring vertices
2677 */
2678 printf("\tVertex neighbors...");
2679 if (s->neighbor_vert) {
2680 for (k = 0; k < s->np; k++)
2681 FREE_17(s->neighbor_vert[k]);
2682 FREE_17(s->neighbor_vert);
2683 }
2684 FREE_17(s->nneighbor_vert);
2685 s->neighbor_vert = MALLOC_17(s->np,int *);
2686 s->nneighbor_vert = MALLOC_17(s->np,int);
2687 /*
2688 * We know the number of neighbors beforehand
2689 */
2690 if (border) {
2691 for (k = 0; k < s->np; k++) {
2692 if (s->nneighbor_tri[k] > 0) {
2693 if (border[k]) {
2694 s->neighbor_vert[k] = MALLOC_17(s->nneighbor_tri[k]+1,int);
2695 s->nneighbor_vert[k] = s->nneighbor_tri[k]+1;
2696 }
2697 else {
2698 s->neighbor_vert[k] = MALLOC_17(s->nneighbor_tri[k],int);
2699 s->nneighbor_vert[k] = s->nneighbor_tri[k];
2700 }
2701 }
2702 else {
2703 s->neighbor_vert[k] = NULL;
2704 s->nneighbor_vert[k] = 0;
2705 }
2706 }
2707 }
2708 else {
2709 for (k = 0; k < s->np; k++) {
2710 if (s->nneighbor_tri[k] > 0) {
2711 s->neighbor_vert[k] = MALLOC_17(s->nneighbor_tri[k],int);
2712 s->nneighbor_vert[k] = s->nneighbor_tri[k];
2713 }
2714 else {
2715 s->neighbor_vert[k] = NULL;
2716 s->nneighbor_vert[k] = 0;
2717 }
2718 }
2719 }
2720 nfix_distinct = 0;
2721 for (k = 0; k < s->np; k++) {
2722 neighbors = s->neighbor_vert[k];
2723 nneighbors = 0;
2724 for (p = 0; p < s->nneighbor_tri[k]; p++) {
2725 /*
2726 * Fit in the other vertices of the neighboring triangle
2727 */
2728 for (c = 0; c < 3; c++) {
2729 vert = s->tris[s->neighbor_tri[k][p]].vert[c];
2730 if (vert != k) {
2731 for (q = 0, found = FALSE; q < nneighbors; q++) {
2732 if (neighbors[q] == vert) {
2733 found = TRUE;
2734 break;
2735 }
2736 }
2737 if (!found) {
2738 if (nneighbors < s->nneighbor_vert[k])
2739 neighbors[nneighbors++] = vert;
2740 else if (!border || !border[k]) {
2741 if (check_too_many_neighbors) {
2742 printf("Too many neighbors for vertex %d.",k);
2743 return FAIL;
2744 }
2745 else
2746 printf("\tWarning: Too many neighbors for vertex %d\n",k);
2747 }
2748 }
2749 }
2750 }
2751 }
2752 if (nneighbors != s->nneighbor_vert[k]) {
2753#ifdef REPORT_WARNINGS
2754 printf("\n\tIncorrect number of distinct neighbors for vertex %d (%d instead of %d) [fixed].",
2755 k,nneighbors,s->nneighbor_vert[k]);
2756#endif
2757 nfix_distinct++;
2758 s->nneighbor_vert[k] = nneighbors;
2759 }
2760 }
2761 printf("[done]\n");
2762 /*
2763 * Distance calculation follows
2764 */
2765 calculate_vertex_distances(s);
2766 mne_compute_surface_cm((MneSurfaceOld*)s);
2767 /*
2768 * Summarize the defects
2769 */
2770 if (nfix_defect > 0)
2771 printf("\tWarning: %d topological defects were fixed.\n",nfix_defect);
2772 if (nfix_distinct > 0)
2773 printf("\tWarning: %d vertices had incorrect number of distinct neighbors (fixed).\n",nfix_distinct);
2774 if (nfix_no_neighbors > 0)
2775 printf("\tWarning: %d vertices did not have any neighboring triangles (fixed)\n",nfix_no_neighbors);
2776#ifdef DEBUG
2777 for (k = 0; k < s->np; k++) {
2778 if (s->nneighbor_vert[k] <= 0)
2779 printf("No neighbors for vertex %d\n",k);
2780 if (s->nneighbor_tri[k] <= 0)
2781 printf("No neighbor tris for vertex %d\n",k);
2782 }
2783#endif
2784 return OK;
2785}
2786
2787//=============================================================================================================
2788
2789int MneSurfaceOrVolume::mne_source_space_add_geometry_info(MneSourceSpaceOld* s, int do_normals)
2790{
2791 return add_geometry_info(s,do_normals,NULL,TRUE);
2792}
2793
2794//=============================================================================================================
2795
2796int MneSurfaceOrVolume::mne_source_space_add_geometry_info2(MneSourceSpaceOld* s, int do_normals)
2797
2798{
2799 return add_geometry_info(s,do_normals,NULL,FALSE);
2800}
2801
2802//=============================================================================================================
2803
2804int MneSurfaceOrVolume::mne_label_area(char *label, MneSourceSpaceOld* s, float *areap) /* Return the area here */
2805/*
2806 * Calculate the area of the label
2807 */
2808{
2809 int *sel = NULL;
2810 int nsel = 0;
2811 float area;
2812 int k,q,nneigh,*neigh;
2813
2814 if (!s) {
2815 qCritical("Source space not specified for mne_label_area");
2816 goto bad;
2817 }
2818 if (mne_read_label(label,NULL,&sel,&nsel))
2819 goto bad;
2820
2821 area = 0.0;
2822 for (k = 0; k < nsel; k++) {
2823 if (sel[k] < 0 || sel[k] >= s->np) {
2824 qCritical("Label vertex index out of range in mne_label_area");
2825 goto bad;
2826 }
2827 nneigh = s->nneighbor_tri[sel[k]];
2828 neigh = s->neighbor_tri[sel[k]];
2829 for (q = 0; q < nneigh; q++)
2830 area += s->tris[neigh[q]].area/3.0;
2831 }
2832 FREE_17(sel);
2833 *areap = area;
2834 return OK;
2835
2836bad : {
2837 FREE_17(sel);
2838 return FAIL;
2839 }
2840}
2841
2842//=============================================================================================================
2843
2844// Align the MEG fiducials to the MRI fiducials
2845int MneSurfaceOrVolume::align_fiducials(FiffDigitizerData* head_dig,
2846 FiffDigitizerData* mri_dig,
2847 MneMshDisplaySurface* head_surf,
2848 int niter,
2849 int scale_head,
2850 float omit_dist,
2851 float *scales)
2852
2853{
2854 float *head_fid[3],*mri_fid[3],**fid;
2855 int j,k;
2856 FiffDigPoint p;
2857 FiffDigitizerData* dig = NULL;
2858 float nasion_weight = 5.0;
2859
2860 if (!head_dig) {
2861 qCritical("MEG head coordinate system digitizer data not available");
2862 goto bad;
2863 }
2864 if (!mri_dig) {
2865 qCritical("MRI coordinate system digitizer data not available");
2866 goto bad;
2867 }
2868
2869 for (j = 0; j < 2; j++) {
2870 dig = j == 0 ? head_dig : mri_dig;
2871 fid = j == 0 ? head_fid : mri_fid;
2872
2873 for (k = 0; k < 3; k++) {
2874 fid[k] = NULL;
2875 }
2876
2877 for (k = 0; k < dig->npoint; k++) {
2878 p = dig->points[k];
2879 if (p.kind == FIFFV_POINT_CARDINAL) {
2880 if (p.ident == FIFFV_POINT_LPA) {
2881 fid[0] = dig->points[k].r;
2882 }
2883 else if (p.ident == FIFFV_POINT_NASION) {
2884 fid[1] = dig->points[k].r;
2885 }
2886 else if (p.ident == FIFFV_POINT_RPA) {
2887 fid[2] = dig->points[k].r;
2888 }
2889 }
2890 }
2891 }
2892
2893 for (k = 0; k < 3; k++) {
2894 if (!head_fid[k]) {
2895 qCritical("Some of the MEG fiducials were missing");
2896 goto bad;
2897 }
2898
2899 if (!mri_fid[k]) {
2900 qCritical("Some of the MRI fiducials were missing");
2901 goto bad;
2902 }
2903 }
2904
2905 if (scale_head) {
2906 get_head_scale(head_dig,mri_fid,head_surf,scales);
2907 printf("xscale = %.3f yscale = %.3f zscale = %.3f\n",scales[0],scales[1],scales[2]);
2908
2909 for (j = 0; j < 3; j++)
2910 for (k = 0; k < 3; k++)
2911 mri_fid[j][k] = mri_fid[j][k]*scales[k];
2912
2913 scale_display_surface(head_surf,scales);
2914 }
2915
2916 // Initial alignment
2917 FREE_17(head_dig->head_mri_t_adj);
2918 head_dig->head_mri_t_adj = FIFFLIB::FiffCoordTransOld::fiff_make_transform_card(FIFFV_COORD_HEAD,FIFFV_COORD_MRI,
2919 mri_fid[0],mri_fid[1],mri_fid[2]);
2920
2921 for (k = 0; k < head_dig->nfids; k++)
2922 VEC_COPY_17(head_dig->mri_fids[k].r,mri_fid[k]);
2923 FiffCoordTransOld::mne_print_coord_transform_label(stdout,QString("After simple alignment : ").toLatin1().data(),head_dig->head_mri_t_adj);
2924
2925 if (omit_dist > 0)
2926 discard_outlier_digitizer_points(head_dig,head_surf,omit_dist);
2927
2928 // Optional iterative refinement
2929 if (niter > 0 && head_surf) {
2930 for (k = 0; k < niter; k++) {
2931 if (iterate_alignment_once(head_dig,head_surf,nasion_weight,mri_fid[1],k == niter-1 && niter > 1) == FAIL)
2932 goto bad;
2933 }
2934
2935 printf("%d / %d iterations done. RMS dist = %7.1f mm\n",k,niter,
2936 1000.0*rms_digitizer_distance(head_dig,head_surf));
2937 FiffCoordTransOld::mne_print_coord_transform_label(stdout,QString("After refinement :").toLatin1().data(),head_dig->head_mri_t_adj);
2938 }
2939
2940 return OK;
2941
2942bad :
2943 return FAIL;
2944}
2945
2946//=============================================================================================================
2947
2948// Simple head size fit
2949void MneSurfaceOrVolume::get_head_scale(FIFFLIB::FiffDigitizerData* dig,
2950 float **mri_fid,
2951 MneMshDisplaySurface* head_surf,
2952 float *scales)
2953{
2954 float **dig_rr = NULL;
2955 float **head_rr = NULL;
2956 int k,ndig,nhead;
2957 float simplex_size = 2e-2;
2958 float r0[3],Rdig,Rscalp;
2959 float LR[3],LN[3],len,norm[3],diff[3];
2960
2961 scales[0] = scales[1] = scales[2] = 1.0;
2962 if (!dig || !head_surf || !mri_fid){
2963 return;
2964 }
2965
2966 dig_rr = MALLOC_17(dig->npoint,float *);
2967 head_rr = MALLOC_17(head_surf->s->np,float *);
2968
2969 // Pick only the points with positive z
2970 for (k = 0, ndig = 0; k < dig->npoint; k++) {
2971 if (dig->points[k].r[Z_17] > 0) {
2972 dig_rr[ndig++] = dig->points[k].r;
2973 }
2974 }
2975
2976 if (!UTILSLIB::Sphere::fit_sphere_to_points(dig_rr,ndig,simplex_size,r0,&Rdig)){
2977 goto out;
2978 }
2979
2980 printf("Polhemus : (%.1f %.1f %.1f) mm R = %.1f mm\n",1000*r0[X_17],1000*r0[Y_17],1000*r0[Z_17],1000*Rdig);
2981
2982 // Pick only the points above the fiducial plane
2983 VEC_DIFF_17(mri_fid[0],mri_fid[2],LR);
2984 VEC_DIFF_17(mri_fid[0],mri_fid[1],LN);
2985 CROSS_PRODUCT_17(LR,LN,norm);
2986 len = VEC_LEN_17(norm);
2987 norm[0] = norm[0]/len;
2988 norm[1] = norm[1]/len;
2989 norm[2] = norm[2]/len;
2990
2991 for (k = 0, nhead = 0; k < head_surf->s->np; k++) {
2992 VEC_DIFF_17(mri_fid[0],head_surf->s->rr[k],diff);
2993 if (VEC_DOT_17(diff,norm) > 0) {
2994 head_rr[nhead++] = head_surf->s->rr[k];
2995 }
2996 }
2997
2998 if (!UTILSLIB::Sphere::fit_sphere_to_points(head_rr,nhead,simplex_size,r0,&Rscalp)) {
2999 goto out;
3000 }
3001
3002 printf("Scalp : (%.1f %.1f %.1f) mm R = %.1f mm\n",1000*r0[X_17],1000*r0[Y_17],1000*r0[Z_17],1000*Rscalp);
3003
3004 scales[0] = scales[1] = scales[2] = Rdig/Rscalp;
3005
3006out : {
3007 FREE_17(dig_rr);
3008 FREE_17(head_rr);
3009 return;
3010 }
3011}
3012
3013//=============================================================================================================
3014
3015int MneSurfaceOrVolume::discard_outlier_digitizer_points(FIFFLIB::FiffDigitizerData* d,
3016 MneMshDisplaySurface* head,
3017 float maxdist)
3018/*
3019 * Discard outlier digitizer points
3020 */
3021{
3022 int discarded = 0;
3023 int k;
3024
3025 if (d && head) {
3026 d->dist_valid = FALSE;
3027 calculate_digitizer_distances(d,head,TRUE,TRUE);
3028 for (k = 0; k < d->npoint; k++) {
3029 d->discard[k] = FALSE;
3030 /*
3031 * Discard unless cardinal landmark or HPI coil
3032 */
3033 if (std::fabs(d->dist[k]) > maxdist &&
3034 d->points[k].kind != FIFFV_POINT_CARDINAL &&
3035 d->points[k].kind != FIFFV_POINT_HPI) {
3036 discarded++;
3037 d->discard[k] = TRUE;
3038 }
3039 }
3040 printf("%d points discarded (maxdist = %6.1f mm).\n",discarded,1000*maxdist);
3041 }
3042 return discarded;
3043}
3044
3045//=============================================================================================================
3046
3047void MneSurfaceOrVolume::calculate_digitizer_distances(FIFFLIB::FiffDigitizerData* dig, MneMshDisplaySurface* head,
3048 int do_all, int do_approx)
3049/*
3050 * Calculate the distances from the scalp surface
3051 */
3052{
3053 float** rr = ALLOC_CMATRIX_17(dig->npoint,3);
3054 int k,nactive;
3055 int* closest;
3056 float* dist;
3057 FiffDigPoint point;
3058 FiffCoordTransOld* t = dig->head_mri_t_adj ? dig->head_mri_t_adj : dig->head_mri_t;
3059 int nstep = 4;
3060
3061 if (dig->dist_valid)
3062 {
3063 FREE_CMATRIX_17(rr);
3064 return ;
3065 }
3066
3067 dig->dist = REALLOC_17(dig->dist,dig->npoint,float);
3068 if (!dig->closest) {
3069 /*
3070 * Ensure that all closest values are initialized correctly
3071 */
3072 dig->closest = MALLOC_17(dig->npoint,int);
3073 for (k = 0; k < dig->npoint; k++)
3074 dig->closest[k] = -1;
3075 }
3076 FREE_CMATRIX_17(dig->closest_point);
3077
3078 dig->closest_point = ALLOC_CMATRIX_17(dig->npoint,3);
3079 closest = MALLOC_17(dig->npoint,int);
3080 dist = MALLOC_17(dig->npoint,float);
3081
3082 for (k = 0, nactive = 0; k < dig->npoint; k++) {
3083 if ((dig->active[k] && !dig->discard[k]) || do_all) {
3084 point = dig->points.at(k);
3085 VEC_COPY_17(rr[nactive],point.r);
3086 FiffCoordTransOld::fiff_coord_trans(rr[nactive],t,FIFFV_MOVE);
3087 if (do_approx) {
3088 closest[nactive] = dig->closest[k];
3089 if (closest[nactive] < 0)
3090 do_approx = FALSE;
3091 }
3092 else
3093 closest[nactive] = -1;
3094 nactive++;
3095 }
3096 }
3097
3098 mne_find_closest_on_surface_approx(head->s,rr,nactive,closest,dist,nstep);
3099 /*
3100 * Project the points on the triangles
3101 */
3102 if (!do_approx)
3103 printf("Inside or outside for %d points...",nactive);
3104 for (k = 0, nactive = 0; k < dig->npoint; k++) {
3105 if ((dig->active[k] && !dig->discard[k]) || do_all) {
3106 dig->dist[k] = dist[nactive];
3107 dig->closest[k] = closest[nactive];
3108 mne_project_to_triangle(head->s,dig->closest[k],rr[nactive],dig->closest_point[k]);
3109 /*
3110 * The above distance is with respect to the closest triangle only
3111 * We need to use the solid angle criterion to decide the sign reliably
3112 */
3113 if (!do_approx && FALSE) {
3114 if (sum_solids(rr[nactive],head->s)/(4*M_PI) > 0.9)
3115 dig->dist[k] = - std::fabs(dig->dist[k]);
3116 else
3117 dig->dist[k] = std::fabs(dig->dist[k]);
3118 }
3119 nactive++;
3120 }
3121 }
3122
3123 if (!do_approx)
3124 printf("[done]\n");
3125
3126 FREE_CMATRIX_17(rr);
3127 FREE_17(closest);
3128 FREE_17(dist);
3129 dig->dist_valid = TRUE;
3130
3131 return;
3132}
3133
3134//=============================================================================================================
3135
3136int MneSurfaceOrVolume::iterate_alignment_once(FIFFLIB::FiffDigitizerData* dig, /* The digitizer data */
3137 MneMshDisplaySurface* head, /* The head surface */
3138 int nasion_weight, /* Weight for the nasion */
3139 float *nasion_mri, /* Fixed correspondence point for the nasion (optional) */
3140 int last_step) /* Is this the last iteration step */
3141/*
3142 * Find the best alignment of the coordinate frames
3143 */
3144{
3145 int res = FAIL;
3146 float **rr_head = NULL;
3147 float **rr_mri = NULL;
3148 float *w = NULL;
3149 int k,nactive;
3150 FiffDigPoint point;
3151 FiffCoordTransOld* t = NULL;
3152 float max_diff = 40e-3;
3153
3154 if (!dig->head_mri_t_adj) {
3155 qCritical()<<"Not adjusting the transformation";
3156 goto out;
3157 }
3158 /*
3159 * Calculate initial distances
3160 */
3161 calculate_digitizer_distances(dig,head,FALSE,TRUE);
3162
3163 /*
3164 * Set up the alignment
3165 */
3166 rr_head = ALLOC_CMATRIX_17(dig->npoint,3);
3167 rr_mri = ALLOC_CMATRIX_17(dig->npoint,3);
3168 w = MALLOC_17(dig->npoint,float);
3169
3170 for (k = 0, nactive = 0; k < dig->npoint; k++) {
3171 if (dig->active[k] && !dig->discard[k]) {
3172 point = dig->points.at(k);
3173 VEC_COPY_17(rr_head[nactive],point.r);
3174 VEC_COPY_17(rr_mri[nactive],dig->closest_point[k]);
3175 /*
3176 * Special handling for the nasion
3177 */
3178 if (point.kind == FIFFV_POINT_CARDINAL &&
3179 point.ident == FIFFV_POINT_NASION) {
3180 w[nactive] = nasion_weight;
3181 if (nasion_mri) {
3182 VEC_COPY_17(rr_mri[nactive],nasion_mri);
3183 VEC_COPY_17(rr_head[nactive],nasion_mri);
3184 FiffCoordTransOld::fiff_coord_trans_inv(rr_head[nactive],
3185 dig->head_mri_t_adj ? dig->head_mri_t_adj : dig->head_mri_t,
3186 FIFFV_MOVE);
3187 }
3188 }
3189 else
3190 w[nactive] = 1.0;
3191 nactive++;
3192 }
3193 }
3194 if (nactive < 3) {
3195 qCritical() << "Not enough points to do the alignment";
3196 goto out;
3197 }
3198 if ((t = FiffCoordTransOld::procrustes_align(FIFFV_COORD_HEAD, FIFFV_COORD_MRI,
3199 rr_head, rr_mri, w, nactive, max_diff)) == NULL)
3200 goto out;
3201
3202 if (dig->head_mri_t_adj)
3203 *dig->head_mri_t_adj = *t;
3204 FREE_17(t);
3205 /*
3206 * Calculate final distances
3207 */
3208 dig->dist_valid = FALSE;
3209 calculate_digitizer_distances(dig,head,FALSE,!last_step);
3210 res = OK;
3211 goto out;
3212
3213out : {
3214 FREE_CMATRIX_17(rr_head);
3215 FREE_CMATRIX_17(rr_mri);
3216 FREE_17(w);
3217 return res;
3218 }
3219}
3220
3221//=============================================================================================================
3222
3223float MneSurfaceOrVolume::rms_digitizer_distance(FIFFLIB::FiffDigitizerData* dig, MneMshDisplaySurface* head)
3224{
3225 float rms;
3226 int k,nactive;
3227
3228 calculate_digitizer_distances(dig,head,FALSE,TRUE);
3229
3230 for (k = 0, rms = 0.0, nactive = 0; k < dig->npoint; k++)
3231 if (dig->active[k] && !dig->discard[k]) {
3232 rms = rms + dig->dist[k]*dig->dist[k];
3233 nactive++;
3234 }
3235 if (nactive > 1)
3236 rms = rms/(nactive-1);
3237 return sqrt(rms);
3238}
3239
3240//=============================================================================================================
3241
3242void MneSurfaceOrVolume::scale_display_surface(MneMshDisplaySurface* surf,
3243 float *scales)
3244/*
3245 * Not quite complete yet
3246 */
3247{
3248 int j,k;
3249
3250 if (!surf || !scales)
3251 return;
3252
3253 for (k = 0; k < 3; k++) {
3254 surf->minv[k] = scales[k]*surf->minv[k];
3255 surf->maxv[k] = scales[k]*surf->maxv[k];
3256 }
3257 for (j = 0; j < surf->s->np; j++)
3258 for (k = 0; k < 3; k++)
3259 surf->s->rr[j][k] = surf->s->rr[j][k]*scales[k];
3260 return;
3261}
3262
3263//=============================================================================================================
3264
3265void MneSurfaceOrVolume::add_uniform_curv(MneSurfaceOld *s)
3266{
3267 int k;
3268 if (!s)
3269 return;
3270 if (s->curv)
3271 return;
3272 s->curv = MALLOC_17(s->np,float);
3273 for (k = 0; k < s->np; k++)
3274 s->curv[k] = 1.0;
3275 return;
3276}
3277
3278//=============================================================================================================
3279
3280char * MneSurfaceOrVolume::mne_compose_surf_name(const char *subj,
3281 const char *name,
3282 const char *prefix)
3283/*
3284 * Get the full path to a surface using the FreeSurfer hierarchy
3285 */
3286{
3287 char *res;
3288 char *subjects_dir = getenv("SUBJECTS_DIR");
3289
3290 if (!subjects_dir || strlen(subjects_dir) == 0) {
3291 qCritical()<<"SUBJECTS_DIR not set. Cannot continue.";
3292 return NULL;
3293 }
3294 if (!subj || strlen(subj) == 0) {
3295 subj = getenv("SUBJECT");
3296 if (subj == NULL || strlen(subj) == 0) {
3297 qCritical()<<"SUBJECT not set. Cannot continue.";
3298 return NULL;
3299 }
3300 }
3301 if (prefix && strlen(prefix) > 0) {
3302 res = MALLOC_17(strlen(subjects_dir)+strlen(subj)+strlen(name)+strlen(prefix)+20,char);
3303 sprintf(res,"%s/%s/surf/%s.%s",subjects_dir,subj,prefix,name);
3304 }
3305 else {
3306 res = MALLOC_17(strlen(subjects_dir)+strlen(subj)+strlen(name)+20,char);
3307 sprintf(res,"%s/%s/surf/%s",subjects_dir,subj,name);
3308 }
3309 return res;
3310}
3311
3312//=============================================================================================================
3313
3314MneSourceSpaceOld* MneSurfaceOrVolume::mne_load_surface(char *surf_file,
3315 char *curv_file)
3316{
3317 return mne_load_surface_geom(surf_file,curv_file,TRUE,TRUE);
3318}
3319
3320//=============================================================================================================
3321
3322MneSourceSpaceOld* MneSurfaceOrVolume::mne_load_surface_geom(char *surf_file,
3323 char *curv_file,
3324 int add_geometry,
3325 int check_too_many_neighbors)
3326 /*
3327 * Load the surface and add the geometry information
3328 */
3329{
3330 float **verts = Q_NULLPTR;
3331 float *curvs = Q_NULLPTR;
3332 int **tris = Q_NULLPTR;
3333 int nvert;
3334 int ntri;
3335 int ncurv;
3336 int k;
3337 MneSourceSpaceOld* s = Q_NULLPTR;
3338 void *tags = Q_NULLPTR;
3339
3340 if (mne_read_triangle_file(surf_file,
3341 &nvert,
3342 &ntri,
3343 &verts,
3344 &tris,
3345 &tags) == -1)
3346 goto bad;
3347
3348 if (curv_file != Q_NULLPTR) {
3349 if (mne_read_curvature_file(curv_file,&curvs,&ncurv) == -1)
3350 goto bad;
3351 if (ncurv != nvert) {
3352 qCritical()<<"Incorrect number of vertices in the curvature file.";
3353 goto bad;
3354 }
3355 }
3356
3357 s = new MneSourceSpaceOld(0);
3358 s->rr = verts; verts = Q_NULLPTR;
3359 s->itris = tris; tris = Q_NULLPTR;
3360 s->ntri = ntri;
3361 s->np = nvert;
3362 s->curv = curvs; curvs = Q_NULLPTR;
3363 s->val = MALLOC_17(s->np,float);
3364 if (add_geometry) {
3365 if (check_too_many_neighbors) {
3366 if (mne_source_space_add_geometry_info(s,TRUE) != OK)
3367 goto bad;
3368 }
3369 else {
3370 if (mne_source_space_add_geometry_info2(s,TRUE) != OK)
3371 goto bad;
3372 }
3373 }
3374 else if (s->nn == Q_NULLPTR) { /* Normals only */
3375 if (mne_add_vertex_normals(s) != OK)
3376 goto bad;
3377 }
3378 else
3379 mne_add_triangle_data(s);
3380 s->nuse = s->np;
3381 s->inuse = MALLOC_17(s->np,int);
3382 s->vertno = MALLOC_17(s->np,int);
3383 for (k = 0; k < s->np; k++) {
3384 s->val[k] = 0.0;
3385 s->inuse[k] = TRUE;
3386 s->vertno[k] = k;
3387 }
3388 s->mgh_tags = tags;
3389 s->vol_geom = mne_get_volume_geom_from_tag(tags);
3390
3391 return s;
3392
3393bad : {
3394 delete ((MneMghTagGroup*)(tags));
3395 FREE_CMATRIX_17(verts);
3396 FREE_17(curvs);
3397 FREE_ICMATRIX_17(tris);
3398 delete s;
3399 return Q_NULLPTR;
3400 }
3401}
3402
3403//=============================================================================================================
3404
3405int MneSurfaceOrVolume::mne_read_triangle_file(char *fname,
3406 int *nvertp,
3407 int *ntrip,
3408 float ***vertp,
3409 int ***trip,
3410 void **tagsp)
3411/*
3412 * Read the FS triangulated surface
3413 */
3414{
3415 FILE *fp = fopen(fname,"r");
3416 int magic;
3417 char c;
3418
3419 qint32 nvert,ntri,nquad;
3420 float **vert = NULL;
3421 int **tri = NULL;
3422 int k,p;
3423 int quad[4];
3424 int val;
3425 float *rr[5];
3426 int which;
3427
3428 if (fp == NULL) {
3429 qCritical("%s", fname);
3430 goto bad;
3431 }
3432 if (mne_read_int3(fp,&magic) != 0) {
3433 printf("Bad magic in %s",fname);
3434 goto bad;
3435 }
3436 if (magic != TRIANGLE_FILE_MAGIC_NUMBER &&
3437 magic != QUAD_FILE_MAGIC_NUMBER &&
3438 magic != NEW_QUAD_FILE_MAGIC_NUMBER) {
3439 printf("Bad magic in %s (%x vs %x)",fname,magic,
3440 TRIANGLE_FILE_MAGIC_NUMBER);
3441 goto bad;
3442 }
3443 if (magic == TRIANGLE_FILE_MAGIC_NUMBER) {
3444 /*
3445 * Get the comment
3446 */
3447 printf("Triangle file : ");
3448 for (c = fgetc(fp); c != '\n'; c = fgetc(fp)) {
3449 if (c == EOF) {
3450 qCritical()<<"Bad triangle file.";
3451 goto bad;
3452 }
3453 putc(c,stderr);
3454 }
3455 c = fgetc(fp);
3456 /*
3457 * How many vertices and triangles?
3458 */
3459 if (mne_read_int(fp,&nvert) != 0)
3460 goto bad;
3461 if (mne_read_int(fp,&ntri) != 0)
3462 goto bad;
3463 printf(" nvert = %d ntri = %d\n",nvert,ntri);
3464 vert = ALLOC_CMATRIX_17(nvert,3);
3465 tri = ALLOC_ICMATRIX_17(ntri,3);
3466 /*
3467 * Read the vertices
3468 */
3469 for (k = 0; k < nvert; k++) {
3470 if (mne_read_float(fp,vert[k]+X_17) != 0)
3471 goto bad;
3472 if (mne_read_float(fp,vert[k]+Y_17) != 0)
3473 goto bad;
3474 if (mne_read_float(fp,vert[k]+Z_17) != 0)
3475 goto bad;
3476 }
3477 /*
3478 * Read the triangles
3479 */
3480 for (k = 0; k < ntri; k++) {
3481 if (mne_read_int(fp,tri[k]+X_17) != 0)
3482 goto bad;
3483 if (check_vertex(tri[k][X_17],nvert) != OK)
3484 goto bad;
3485 if (mne_read_int(fp,tri[k]+Y_17) != 0)
3486 goto bad;
3487 if (check_vertex(tri[k][Y_17],nvert) != OK)
3488 goto bad;
3489 if (mne_read_int(fp,tri[k]+Z_17) != 0)
3490 goto bad;
3491 if (check_vertex(tri[k][Z_17],nvert) != OK)
3492 goto bad;
3493 }
3494 }
3495 else if (magic == QUAD_FILE_MAGIC_NUMBER ||
3496 magic == NEW_QUAD_FILE_MAGIC_NUMBER) {
3497 if (mne_read_int3(fp,&nvert) != 0)
3498 goto bad;
3499 if (mne_read_int3(fp,&nquad) != 0)
3500 goto bad;
3501 printf("%s file : nvert = %d nquad = %d\n",
3502 magic == QUAD_FILE_MAGIC_NUMBER ? "Quad" : "New quad",
3503 nvert,nquad);
3504 vert = ALLOC_CMATRIX_17(nvert,3);
3505 if (magic == QUAD_FILE_MAGIC_NUMBER) {
3506 for (k = 0; k < nvert; k++) {
3507 if (mne_read_int2(fp,&val) != 0)
3508 goto bad;
3509 vert[k][X_17] = val/100.0;
3510 if (mne_read_int2(fp,&val) != 0)
3511 goto bad;
3512 vert[k][Y_17] = val/100.0;
3513 if (mne_read_int2(fp,&val) != 0)
3514 goto bad;
3515 vert[k][Z_17] = val/100.0;
3516 }
3517 }
3518 else { /* NEW_QUAD_FILE_MAGIC_NUMBER */
3519 for (k = 0; k < nvert; k++) {
3520 if (mne_read_float(fp,vert[k]+X_17) != 0)
3521 goto bad;
3522 if (mne_read_float(fp,vert[k]+Y_17) != 0)
3523 goto bad;
3524 if (mne_read_float(fp,vert[k]+Z_17) != 0)
3525 goto bad;
3526 }
3527 }
3528 ntri = 2*nquad;
3529 tri = ALLOC_ICMATRIX_17(ntri,3);
3530 for (k = 0, ntri = 0; k < nquad; k++) {
3531 for (p = 0; p < 4; p++) {
3532 if (mne_read_int3(fp,quad+p) != 0)
3533 goto bad;
3534 rr[p] = vert[quad[p]];
3535 }
3536 rr[4] = vert[quad[0]];
3537 if (check_quad(rr) != OK)
3538 goto bad;
3539
3540 /*
3541 * The randomization is borrowed from FreeSurfer code
3542 * Strange...
3543 */
3544#define EVEN(n) ((((n) / 2) * 2) == n)
3545#ifdef FOO
3546#define WHICH_FACE_SPLIT(vno0, vno1) \
3547 (1*nearbyint(sqrt(1.9*vno0) + sqrt(3.5*vno1)))
3548
3549 which = WHICH_FACE_SPLIT(quad[0], quad[1]) ;
3550#endif
3551 which = quad[0];
3552 /*
3553 printf("%f ",sqrt(1.9*quad[0]) + sqrt(3.5*quad[1]));
3554 */
3555
3556 if (EVEN(which)) {
3557 tri[ntri][X_17] = quad[0];
3558 tri[ntri][Y_17] = quad[1];
3559 tri[ntri][Z_17] = quad[3];
3560 ntri++;
3561
3562 tri[ntri][X_17] = quad[2];
3563 tri[ntri][Y_17] = quad[3];
3564 tri[ntri][Z_17] = quad[1];
3565 ntri++;
3566 }
3567 else {
3568 tri[ntri][X_17] = quad[0];
3569 tri[ntri][Y_17] = quad[1];
3570 tri[ntri][Z_17] = quad[2];
3571 ntri++;
3572
3573 tri[ntri][X_17] = quad[0];
3574 tri[ntri][Y_17] = quad[2];
3575 tri[ntri][Z_17] = quad[3];
3576 ntri++;
3577 }
3578 }
3579 }
3580 /*
3581 * Optionally read the tags
3582 */
3583 if (tagsp) {
3584 void *tags = NULL;
3585 if (mne_read_mgh_tags(fp, &tags) == FAIL) {
3586 delete((MneMghTagGroup*)tags);
3587 goto bad;
3588 }
3589 *tagsp = tags;
3590 }
3591 fclose(fp);
3592 *nvertp = nvert;
3593 *ntrip = ntri;
3594 *vertp = vert;
3595 *trip = tri;
3596 for (k = 0; k < nvert; k++) {
3597 vert[k][X_17] = vert[k][X_17]/1000.0;
3598 vert[k][Y_17] = vert[k][Y_17]/1000.0;
3599 vert[k][Z_17] = vert[k][Z_17]/1000.0;
3600 }
3601#ifdef FOO
3602 /*
3603 * Which ordering does the file have???
3604 */
3605 for (k = 0; k < ntri; k++) {
3606 help = tri[k][X_17];
3607 tri[k][X_17] = tri[k][Y_17];
3608 tri[k][Y_17] = help;
3609 }
3610#endif
3611 return OK;
3612
3613bad : {
3614 if (fp)
3615 fclose(fp);
3616 FREE_CMATRIX_17(vert);
3617 FREE_ICMATRIX_17(tri);
3618 return FAIL;
3619 }
3620}
3621
3622//=============================================================================================================
3623
3624int MneSurfaceOrVolume::mne_read_curvature_file(char *fname,
3625 float **curvsp,
3626 int *ncurvp)
3627
3628{
3629 FILE *fp = fopen(fname,"r");
3630 int magic;
3631
3632 float *curvs = NULL;
3633 float curvmin,curvmax;
3634 int ncurv = 0;
3635 int nface,val_pervert;
3636 int val,k;
3637
3638 if (!fp) {
3639 qCritical("%s", fname);
3640 goto bad;
3641 }
3642 if (mne_read_int3(fp,&magic) != 0) {
3643 printf( "Bad magic in %s",fname);
3644 goto bad;
3645 }
3646 if (magic == CURVATURE_FILE_MAGIC_NUMBER) { /* A new-style curvature file */
3647 /*
3648 * How many and faces
3649 */
3650 if (mne_read_int(fp,&ncurv) != 0)
3651 goto bad;
3652 if (mne_read_int(fp,&nface) != 0)
3653 goto bad;
3654#ifdef DEBUG
3655 printf("nvert = %d nface = %d\n",ncurv,nface);
3656#endif
3657 if (mne_read_int(fp,&val_pervert) != 0)
3658 goto bad;
3659 if (val_pervert != 1) {
3660 qCritical("Values per vertex not equal to one.");
3661 goto bad;
3662 }
3663 /*
3664 * Read the curvature values
3665 */
3666 curvs = MALLOC_17(ncurv,float);
3667 curvmin = curvmax = 0.0;
3668 for (k = 0; k < ncurv; k++) {
3669 if (mne_read_float(fp,curvs+k) != 0)
3670 goto bad;
3671 if (curvs[k] > curvmax)
3672 curvmax = curvs[k];
3673 if (curvs[k] < curvmin)
3674 curvmin = curvs[k];
3675 }
3676 }
3677 else { /* An old-style curvature file */
3678 ncurv = magic;
3679 /*
3680 * How many vertices
3681 */
3682 if (mne_read_int3(fp,&nface) != 0)
3683 goto bad;
3684#ifdef DEBUG
3685 printf("nvert = %d nface = %d\n",ncurv,nface);
3686#endif
3687 /*
3688 * Read the curvature values
3689 */
3690 curvs = MALLOC_17(ncurv,float);
3691 curvmin = curvmax = 0.0;
3692 for (k = 0; k < ncurv; k++) {
3693 if (mne_read_int2(fp,&val) != 0)
3694 goto bad;
3695 curvs[k] = (float)val/100.0;
3696 if (curvs[k] > curvmax)
3697 curvmax = curvs[k];
3698 if (curvs[k] < curvmin)
3699 curvmin = curvs[k];
3700
3701 }
3702 }
3703#ifdef DEBUG
3704 printf("Curvature range: %f...%f\n",curvmin,curvmax);
3705#endif
3706 *ncurvp = ncurv;
3707 *curvsp = curvs;
3708 return OK;
3709
3710bad : {
3711 if (fp)
3712 fclose(fp);
3713 FREE_17(curvs);
3714 return FAIL;
3715 }
3716}
3717
3718//=============================================================================================================
3719
3720int MneSurfaceOrVolume::check_quad(float **rr)
3721
3722{
3723 float diff[3];
3724 float size;
3725 int k;
3726
3727 return OK;
3728
3729 for (k = 0; k < 4; k++) {
3730 VEC_DIFF_17(rr[k],rr[k+1],diff);
3731 size = VEC_LEN_17(diff);
3732 if (size < 0.1) {
3733 printf("Degenerate quad found. size length = %f mm",size);
3734 return FAIL;
3735 }
3736 }
3737 return OK;
3738}
3739
3740//=============================================================================================================
3741
3742int MneSurfaceOrVolume::check_vertex(int no, int maxno)
3743
3744{
3745 if (no < 0 || no > maxno-1) {
3746 printf("Illegal vertex number %d (max %d).",no,maxno);
3747 return FAIL;
3748 }
3749 return OK;
3750}
3751
3752//=============================================================================================================
3753
3754MneVolGeom* MneSurfaceOrVolume::mne_get_volume_geom_from_tag(void *tagsp)
3755{
3756 MneMghTagGroup* tags = (MneMghTagGroup*)tagsp;
3757 MneMghTag* tag = NULL;
3758 MneVolGeom* vg = NULL;
3759 int k;
3760
3761 if (tags) {
3762 for (k = 0; k < tags->ntags; k++)
3763 if (tags->tags[k]->tag == TAG_OLD_SURF_GEOM) {
3764 tag = tags->tags[k];
3765 break;
3766 }
3767 if (tag)
3768 vg = mne_dup_vol_geom((MneVolGeom*)tag->data);
3769 }
3770 return vg;
3771}
3772
3773//=============================================================================================================
3774
3775MneVolGeom* MneSurfaceOrVolume::mne_dup_vol_geom(MneVolGeom* g)
3776{
3777 MneVolGeom* dup = NULL;
3778 if (g) {
3779 dup = new MneVolGeom();
3780 *dup = *g;
3781 dup->filename = g->filename;
3782 }
3783 return dup;
3784}
3785
3786//=============================================================================================================
3787
3788int MneSurfaceOrVolume::mne_read_mgh_tags(FILE *fp, void **tagsp)
3789/*
3790 * Read all the tags from the file
3791 */
3792{
3793 long long len;
3794 int tag;
3795 unsigned char *tag_data;
3796 MneMghTagGroup **tags = (MneMghTagGroup **)tagsp;
3797
3798 while (1) {
3799 if (read_next_tag(fp,&tag,&len,&tag_data) == FAIL)
3800 return FAIL;
3801 if (tag == 0)
3802 break;
3803 *tags = mne_add_mgh_tag_to_group(*tags,tag,len,tag_data);
3804 }
3805 tagsp = (void **)tags;
3806 return OK;
3807}
3808
3809//=============================================================================================================
3810
3811int MneSurfaceOrVolume::read_next_tag(FILE *fp, int *tagp, long long *lenp, unsigned char **datap)
3812/*
3813 * Read the next tag in the file
3814 */
3815{
3816 int ilen,tag;
3817 long long len;
3818
3819 if (mne_read_int(fp,&tag) == FAIL) {
3820 *tagp = 0;
3821 return OK;
3822 }
3823 if (feof(fp)) {
3824 *tagp = 0;
3825 return OK;
3826 }
3827 switch (tag) {
3828 case TAG_OLD_MGH_XFORM: /* This is obviously a burden of the past */
3829 if (mne_read_int(fp,&ilen) == FAIL)
3830 return FAIL;
3831 len = ilen - 1;
3832 break ;
3833 case TAG_OLD_SURF_GEOM:
3834 case TAG_OLD_USEREALRAS:
3835 case TAG_OLD_COLORTABLE:
3836 len = 0 ;
3837 break ;
3838 default:
3839 if (mne_read_long(fp,&len) == FAIL)
3840 return FAIL;
3841 break;
3842 }
3843 *lenp = len;
3844 *tagp = tag;
3845 if (read_tag_data(fp,tag,len,datap,lenp) == FAIL)
3846 return FAIL;
3847 return OK;
3848}
3849
3850//=============================================================================================================
3851
3852int MneSurfaceOrVolume::read_tag_data(FILE *fp, int tag, long long nbytes, unsigned char **val, long long *nbytesp)
3853/*
3854 * Read the data of one tag
3855 */
3856{
3857 unsigned char *dum = NULL;
3858 size_t snbytes = nbytes;
3859
3860 *val = NULL;
3861 if (nbytes > 0) {
3862 dum = MALLOC_17(nbytes+1,unsigned char);
3863 if (fread(dum,sizeof(unsigned char),nbytes,fp) != snbytes) {
3864 printf("Failed to read %d bytes of tag data",static_cast<int>(nbytes));
3865 FREE_17(dum);
3866 return FAIL;
3867 }
3868 dum[nbytes] = '\0'; /* Ensure null termination */
3869 *val = dum;
3870 *nbytesp = nbytes;
3871 }
3872 else { /* Need to handle special cases */
3873 if (tag == TAG_OLD_SURF_GEOM) {
3874 MneVolGeom* g = read_vol_geom(fp);
3875 if (!g)
3876 return FAIL;
3877 *val = (unsigned char *)g;
3878 *nbytesp = sizeof(MneVolGeom);
3879 }
3880 else if (tag == TAG_OLD_USEREALRAS || tag == TAG_USEREALRAS) {
3881 int *vi = MALLOC_17(1,int);
3882 if (mne_read_int(fp,vi) == FAIL)
3883 vi = 0;
3884 *val = (unsigned char *)vi;
3885 *nbytesp = sizeof(int);
3886 }
3887 else {
3888 printf("Encountered an unknown tag with no length specification : %d\n",tag);
3889 *val = NULL;
3890 *nbytesp = 0;
3891 }
3892 }
3893 return OK;
3894}
3895
3896//=============================================================================================================
3897
3898MneMghTagGroup* MneSurfaceOrVolume::mne_add_mgh_tag_to_group(MneMghTagGroup* g, int tag, long long len, unsigned char *data)
3899{
3900 MneMghTag* new_tag = NULL;
3901
3902 if (!g)
3903 g = new MneMghTagGroup();
3904 g->tags = REALLOC_17(g->tags,g->ntags+1,MneMghTag*);
3905 g->tags[g->ntags++] = new_tag = new MneMghTag();
3906 new_tag->tag = tag;
3907 new_tag->len = len;
3908 new_tag->data = data;
3909
3910 return g;
3911}
3912
3913//=============================================================================================================
3914
3915MneVolGeom* MneSurfaceOrVolume::read_vol_geom(FILE *fp)
3916/*
3917 * This the volume geometry reading code from FreeSurfer
3918 */
3919{
3920 char line[256];
3921 char param[64];
3922 char eq[2];
3923 char buf[256];
3924 int vgRead = 0;
3925 char *p = 0;
3926 int counter = 0;
3927 long pos = 0;
3928 int fail = 0;
3929
3930 MneVolGeom* vg = new MneVolGeom();
3931
3932 while ((p = fgets(line, sizeof(line), fp)) && counter < 8)
3933 {
3934 if (strlen(p) == 0)
3935 break ;
3936 sscanf(line, "%s %s %*s", param, eq);
3937 if (!strcmp(param, "valid")) {
3938 sscanf(line, "%s %s %d \n", param, eq, &vg->valid);
3939 vgRead = 1;
3940 counter++;
3941 }
3942 else if (!strcmp(param, "filename")) {
3943 if (sscanf(line, "%s %s %s\n", param, eq, buf) >= 3)
3944 vg->filename = mne_strdup(buf);
3945 counter++;
3946 }
3947 else if (!strcmp(param, "volume")) {
3948 sscanf(line, "%s %s %d %d %d\n",
3949 param, eq, &vg->width, &vg->height, &vg->depth);
3950 counter++;
3951 }
3952 else if (!strcmp(param, "voxelsize")) {
3953 sscanf(line, "%s %s %f %f %f\n",
3954 param, eq, &vg->xsize, &vg->ysize, &vg->zsize);
3955 /*
3956 * We like these to be in meters
3957 */
3958 vg->xsize = vg->xsize/1000.0;
3959 vg->ysize = vg->ysize/1000.0;
3960 vg->zsize = vg->zsize/1000.0;
3961 counter++;
3962 }
3963 else if (!strcmp(param, "xras")) {
3964 sscanf(line, "%s %s %f %f %f\n",
3965 param, eq, vg->x_ras, vg->x_ras+1, vg->x_ras+2);
3966 counter++;
3967 }
3968 else if (!strcmp(param, "yras")) {
3969 sscanf(line, "%s %s %f %f %f\n",
3970 param, eq, vg->y_ras, vg->y_ras+1, vg->y_ras+2);
3971 counter++;
3972 }
3973 else if (!strcmp(param, "zras")) {
3974 sscanf(line, "%s %s %f %f %f\n",
3975 param, eq, vg->z_ras, vg->z_ras+1, vg->z_ras+2);
3976 counter++;
3977 }
3978 else if (!strcmp(param, "cras")) {
3979 sscanf(line, "%s %s %f %f %f\n",
3980 param, eq, vg->c_ras, vg->c_ras+1, vg->c_ras+2);
3981 vg->c_ras[0] = vg->c_ras[0]/1000.0;
3982 vg->c_ras[1] = vg->c_ras[1]/1000.0;
3983 vg->c_ras[2] = vg->c_ras[2]/1000.0;
3984 counter++;
3985 }
3986 /* rememver the current position */
3987 pos = ftell(fp); /* if fail = 0, then ok */
3988 };
3989 if (p) { /* we read one more line */
3990 if (pos > 0 ) /* if success in getting pos, then */
3991 fail = fseek(fp, pos, SEEK_SET); /* restore the position */
3992 /* note that this won't allow compression using pipe */
3993 }
3994 if (!vgRead) {
3995 delete vg;
3996 vg = new MneVolGeom();
3997 }
3998 (void)fail; // squash compiler warning, set but unused
3999 return vg;
4000}
4001
4002//=============================================================================================================
4003
4004int MneSurfaceOrVolume::mne_read_int3(FILE *in, int *ival)
4005/*
4006 * Read the strange 3-byte integer
4007 */
4008{
4009 unsigned int s = 0;
4010
4011 if (fread (&s,3,1,in) != 1) {
4012 if (ferror(in))
4013 qCritical("mne_read_int3");
4014 else
4015 qCritical("mne_read_int3 could not read data");
4016 return FAIL;
4017 }
4018 s = (unsigned int)UTILSLIB::IOUtils::swap_int(s);
4019 *ival = ((s >> 8) & 0xffffff);
4020 return OK;
4021}
4022
4023//=============================================================================================================
4024
4025int MneSurfaceOrVolume::mne_read_int(FILE *in, qint32 *ival)
4026/*
4027 * Read a 32-bit integer
4028 */
4029{
4030 qint32 s ;
4031 if (fread (&s,sizeof(qint32),1,in) != 1) {
4032 if (ferror(in))
4033 qCritical("mne_read_int");
4034 else
4035 qCritical("mne_read_int could not read data");
4036 return FAIL;
4037 }
4038 *ival = UTILSLIB::IOUtils::swap_int(s);
4039 return OK;
4040}
4041
4042//=============================================================================================================
4043
4044int MneSurfaceOrVolume::mne_read_int2(FILE *in, int *ival)
4045/*
4046 * Read int from short
4047 */
4048{
4049 short s ;
4050 if (fread (&s,sizeof(short),1,in) != 1) {
4051 if (ferror(in))
4052 qCritical("mne_read_int2");
4053 else
4054 qCritical("mne_read_int2 could not read data");
4055 return FAIL;
4056 }
4057 *ival = UTILSLIB::IOUtils::swap_short(s);
4058 return OK;
4059}
4060
4061//=============================================================================================================
4062
4063int MneSurfaceOrVolume::mne_read_float(FILE *in, float *fval)
4064/*
4065 * Read float
4066 */
4067{
4068 float f ;
4069 if (fread (&f,sizeof(float),1,in) != 1) {
4070 if (ferror(in))
4071 qCritical("mne_read_float");
4072 else
4073 qCritical("mne_read_float could not read data");
4074 return FAIL;
4075 }
4076 *fval = UTILSLIB::IOUtils::swap_float(f);
4077 return OK;
4078}
4079
4080//=============================================================================================================
4081
4082int MneSurfaceOrVolume::mne_read_long(FILE *in, long long *lval)
4083/*
4084 * Read a 64-bit integer
4085 */
4086{
4087 long long s ;
4088 if (fread (&s,sizeof(long long),1,in) != 1) {
4089 if (ferror(in))
4090 qCritical("mne_read_long");
4091 else
4092 qCritical("mne_read_long could not read data");
4093 return FAIL;
4094 }
4095 *lval = UTILSLIB::IOUtils::swap_long(s);
4096 return OK;
4097}
4098
4099//=============================================================================================================
4100
4101char *MneSurfaceOrVolume::mne_strdup(const char *s)
4102{
4103 char *res;
4104 if (s == NULL)
4105 return NULL;
4106 res = MALLOC_17(strlen(s)+1,char);
4107 strcpy(res,s);
4108 return res;
4109}
4110
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_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
FIFFV_MNE_COORD_MRI_VOXEL
#define FIFFV_MNE_COORD_MRI_VOXEL
Definition fiff_constants.h:482
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
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_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
FIFFV_MNE_COORD_RAS
#define FIFFV_MNE_COORD_RAS
Definition fiff_constants.h:483
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_stream.h
FiffStream class declaration.
fiff_digitizer_data.h
FiffDigitizerData class declaration.
k
int k
Definition fiff_tag.cpp:324
FIFFB_BEM_SURF
#define FIFFB_BEM_SURF
Definition fiff_file.h:404
FIFF_BEM_SURF_NNODE
#define FIFF_BEM_SURF_NNODE
Definition fiff_file.h:733
FIFF_BEM_SURF_NODES
#define FIFF_BEM_SURF_NODES
Definition fiff_file.h:735
FIFF_BEM_SURF_TRIANGLES
#define FIFF_BEM_SURF_TRIANGLES
Definition fiff_file.h:736
FIFF_BEM_COORD_FRAME
#define FIFF_BEM_COORD_FRAME
Definition fiff_file.h:743
FIFFB_BEM
#define FIFFB_BEM
Definition fiff_file.h:403
FIFF_BEM_SURF_ID
#define FIFF_BEM_SURF_ID
Definition fiff_file.h:731
FIFF_BEM_SURF_NTRI
#define FIFF_BEM_SURF_NTRI
Definition fiff_file.h:734
FIFF_BEM_SIGMA
#define FIFF_BEM_SIGMA
Definition fiff_file.h:744
FIFF_BEM_SURF_NORMALS
#define FIFF_BEM_SURF_NORMALS
Definition fiff_file.h:737
FIFF_SUBJ_HIS_ID
#define FIFF_SUBJ_HIS_ID
Definition fiff_file.h:573
fiff_dig_point.h
FiffDigPoint class declaration.
ioutils.h
IOUtils class declaration.
sphere.h
Sphere class declaration.
mne_triangle.h
MneTriangle class declaration.
mne_mgh_tag_group.h
MneMghTagGroup class declaration.
mne_vol_geom.h
MneVolGeom class declaration.
filter_thread_arg.h
Filter Thread Argument (FilterThreadArg) class declaration.
mne_patch_info.h
MNE Patch Information (MnePatchInfo) class declaration.
mne_mgh_tag.h
MneMghTag class declaration.
mne_source_space_old.h
MneSourceSpaceOld class declaration.
mne_msh_display_surface.h
MneMshDisplaySurface class declaration.
mne_proj_data.h
MneProjData class declaration.
mne_nearest.h
MneNearest class declaration.
mne_surface_old.h
MneSurfaceOld class declaration.
mne_surface_or_volume.h
MNE Surface or Volume (MneSurfaceOrVolume) class declaration.
FIFFLIB::FiffCoordTransOld
Coordinate transformation descriptor.
Definition fiff_coord_trans_old.h:81
FIFFLIB::FiffDigitizerData
Digitization points container and description.
Definition fiff_digitizer_data.h:74
FIFFLIB::FiffSparseMatrix
Data associated with MNE computations for each mneMeasDataSet.
Definition fiff_sparse_matrix.h:75
FIFFLIB::FiffDigPoint
Digitization point description.
Definition fiff_dig_point.h:69
FIFFLIB::FiffDirNode::SPtr
QSharedPointer< FiffDirNode > SPtr
Definition fiff_dir_node.h:76
FIFFLIB::FiffStream
FIFF File I/O routines.
Definition fiff_stream.h:105
FIFFLIB::FiffStream::SPtr
QSharedPointer< FiffStream > SPtr
Definition fiff_stream.h:107
FIFFLIB::FiffTag::SPtr
QSharedPointer< FiffTag > SPtr
Definition fiff_tag.h:152
MNELIB::FilterThreadArg
Filter Thread Argument Description.
Definition filter_thread_arg.h:77
MNELIB::MneMghTag
The MneMghTag class.
Definition mne_mgh_tag.h:77
MNELIB::MneMghTagGroup
The MneMghTagGroup class.
Definition mne_mgh_tag_group.h:79
MNELIB::MneMshDisplaySurface
The MNE Msh Display Surface class holds information about a surface to be rendered.
Definition mne_msh_display_surface.h:80
MNELIB::MneNearest
This is used in the patch definitions.
Definition mne_nearest.h:78
MNELIB::MnePatchInfo
One item in a derivation data set.
Definition mne_patch_info.h:78
MNELIB::MnePatchInfo::calculate_patch_area
static void calculate_patch_area(MneSourceSpaceOld *s, MnePatchInfo *p)
Definition mne_patch_info.cpp:82
MNELIB::MnePatchInfo::calculate_normal_stats
static void calculate_normal_stats(MneSourceSpaceOld *s, MnePatchInfo *p)
Definition mne_patch_info.cpp:99
MNELIB::MneProjData
The MneProjData class.
Definition mne_proj_data.h:75
MNELIB::MneSourceSpaceOld
This defines a source space.
Definition mne_source_space_old.h:77
MNELIB::MneSurfaceOld
This defines a surface.
Definition mne_surface_old.h:77
MNELIB::MneTriangle
Triangle data.
Definition mne_triangle.h:76
MNELIB::MneVolGeom
MRI data volume geometry information like FreeSurfer keeps it.
Definition mne_vol_geom.h:77
UTILSLIB::IOUtils::swap_int
static qint32 swap_int(qint32 source)
Definition ioutils.cpp:128
UTILSLIB::IOUtils::swap_long
static qint64 swap_long(qint64 source)
Definition ioutils.cpp:162
UTILSLIB::IOUtils::swap_float
static float swap_float(float source)
Definition ioutils.cpp:207
UTILSLIB::IOUtils::swap_short
static qint16 swap_short(qint16 source)
Definition ioutils.cpp:115
UTILSLIB::Sphere::fit_sphere_to_points
static bool fit_sphere_to_points(const Eigen::MatrixXf &rr, float simplex_size, Eigen::VectorXf &r0, float &R)
Generated on Sat Feb 14 2026 16:30:24 for MNE-CPP by doxygen 1.9.8