MNE-CPP  0.1.9
A Framework for Electrophysiology
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"
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>
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 
83 using namespace Eigen;
84 using namespace FIFFLIB;
85 using 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 
139 static 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 
156 float** 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 
187 void mne_free_cmatrix_17 (float **m)
188 {
189  if (m) {
190  FREE_17(*m);
191  FREE_17(m);
192  }
193 }
194 
195 void 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 
219 int **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
235 Eigen::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 
246 void 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 
253 void 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 
258 void 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 
265 void 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 
272 static 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 {
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 
300 static 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 
314 static 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 
328 void 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 
340 MneSurfaceOrVolume::MneSurfaceOrVolume()
341 {
342 }
343 
344 //=============================================================================================================
345 
346 MneSurfaceOrVolume::~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 
403 double 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 
430 double 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 
443 int 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 {
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  if (omit_outside > 0)
523  printf("%d source space points omitted because they are outside the inner skull surface.\n",
524  omit_outside);
525  if (omit > 0)
526  printf("%d source space points omitted because of the %6.1f-mm distance limit.\n",
527  omit,1000*limit);
528  printf("Thank you for waiting.\n");
529  return OK;
530 }
531 
532 //=============================================================================================================
533 
534 int MneSurfaceOrVolume::mne_add_patch_stats(MneSourceSpaceOld* s)
535 {
536  MneNearest* nearest = s->nearest;
537  MneNearest* this_patch;
538  MnePatchInfo* *pinfo = MALLOC_17(s->nuse,MnePatchInfo*);
539  int nave,p,q,k;
540 
541  printf("Computing patch statistics...\n");
542  if (!s->neighbor_tri)
543  if (mne_source_space_add_geometry_info(s,FALSE) != OK)
544  goto bad;
545 
546  if (s->nearest == NULL) {
547  qCritical("The patch information is not available.");
548  goto bad;
549  }
550  if (s->nuse == 0) {
551  FREE_17(s->patches);
552  s->patches = NULL;
553  s->npatch = 0;
554  return OK;
555  }
556  /*
557  * Calculate the average normals and the patch areas
558  */
559  printf("\tareas, average normals, and mean deviations...");
560  mne_sort_nearest_by_nearest(nearest,s->np);
561  nave = 1;
562  for (p = 1, q = 0; p < s->np; p++) {
563  if (nearest[p].nearest != nearest[p-1].nearest) {
564  if (nave == 0) {
565  qCritical("No vertices belong to the patch of vertex %d",nearest[p-1].nearest);
566  goto bad;
567  }
568  if (s->vertno[q] == nearest[p-1].nearest) { /* Some source space points may have been omitted since
569  * the patch information was computed */
570  pinfo[q] = new MnePatchInfo();
571  pinfo[q]->vert = nearest[p-1].nearest;
572  this_patch = nearest+p-nave;
573  pinfo[q]->memb_vert = MALLOC_17(nave,int);
574  pinfo[q]->nmemb = nave;
575  for (k = 0; k < nave; k++) {
576  pinfo[q]->memb_vert[k] = this_patch[k].vert;
577  this_patch[k].patch = pinfo[q];
578  }
579  MnePatchInfo::calculate_patch_area(s,pinfo[q]);
580  MnePatchInfo::calculate_normal_stats(s,pinfo[q]);
581  q++;
582  }
583  nave = 0;
584  }
585  nave++;
586  }
587  if (nave == 0) {
588  qCritical("No vertices belong to the patch of vertex %d",nearest[p-1].nearest);
589  goto bad;
590  }
591  if (s->vertno[q] == nearest[p-1].nearest) {
592  pinfo[q] = new MnePatchInfo;
593  pinfo[q]->vert = nearest[p-1].nearest;
594  this_patch = nearest+p-nave;
595  pinfo[q]->memb_vert = MALLOC_17(nave,int);
596  pinfo[q]->nmemb = nave;
597  for (k = 0; k < nave; k++) {
598  pinfo[q]->memb_vert[k] = this_patch[k].vert;
599  this_patch[k].patch = pinfo[q];
600  }
601  MnePatchInfo::calculate_patch_area(s,pinfo[q]);
602  MnePatchInfo::calculate_normal_stats(s,pinfo[q]);
603  q++;
604  }
605  printf(" %d/%d [done]\n",q,s->nuse);
606 
607  if (s->patches) {
608  for (k = 0; k < s->npatch; k++)
609  if(s->patches[k])
610  delete s->patches[k];
611  FREE_17(s->patches);
612  }
613  s->patches = pinfo;
614  s->npatch = s->nuse;
615 
616  return OK;
617 
618 bad : {
619  FREE_17(pinfo);
620  return FAIL;
621  }
622 }
623 
624 //=============================================================================================================
625 
626 void MneSurfaceOrVolume::rearrange_source_space(MneSourceSpaceOld* s)
627 {
628  int k,p;
629 
630  for (k = 0, s->nuse = 0; k < s->np; k++)
631  if (s->inuse[k])
632  s->nuse++;
633 
634  if (s->nuse == 0) {
635  FREE_17(s->vertno);
636  s->vertno = NULL;
637  }
638  else {
639  s->vertno = REALLOC_17(s->vertno,s->nuse,int);
640  for (k = 0, p = 0; k < s->np; k++)
641  if (s->inuse[k])
642  s->vertno[p++] = k;
643  }
644  if (s->nearest)
645  mne_add_patch_stats(s);
646  return;
647 }
648 
649 //=============================================================================================================
650 
651 void *MneSurfaceOrVolume::filter_source_space(void *arg)
652 {
653  FilterThreadArg* a = (FilterThreadArg*)arg;
654  int p1,p2;
655  double tot_angle;
656  int omit,omit_outside;
657  float r1[3];
658  float mindist,dist,diff[3];
659  int minnode;
660 
661  omit = 0;
662  omit_outside = 0;
663 
664  for (p1 = 0; p1 < a->s->np; p1++) {
665  if (a->s->inuse[p1]) {
666  VEC_COPY_17(r1,a->s->rr[p1]); /* Transform the point to MRI coordinates */
667  if (a->s->coord_frame == FIFFV_COORD_HEAD)
668  FiffCoordTransOld::fiff_coord_trans_inv(r1,a->mri_head_t,FIFFV_MOVE);
669  /*
670  * Check that the source is inside the inner skull surface
671  */
672  tot_angle = sum_solids(r1,a->surf)/(4*M_PI);
673  if (std::fabs(tot_angle-1.0) > 1e-5) {
674  omit_outside++;
675  a->s->inuse[p1] = FALSE;
676  a->s->nuse--;
677  if (a->filtered)
678  fprintf(a->filtered,"%10.3f %10.3f %10.3f\n",
679  1000*r1[X_17],1000*r1[Y_17],1000*r1[Z_17]);
680  }
681  else if (a->limit > 0.0) {
682  /*
683  * Check the distance limit
684  */
685  mindist = 1.0;
686  minnode = 0;
687  for (p2 = 0; p2 < a->surf->np; p2++) {
688  VEC_DIFF_17(r1,a->surf->rr[p2],diff);
689  dist = VEC_LEN_17(diff);
690  if (dist < mindist) {
691  mindist = dist;
692  minnode = p2;
693  }
694  }
695  if (mindist < a->limit) {
696  omit++;
697  a->s->inuse[p1] = FALSE;
698  a->s->nuse--;
699  if (a->filtered)
700  fprintf(a->filtered,"%10.3f %10.3f %10.3f\n",
701  1000*r1[X_17],1000*r1[Y_17],1000*r1[Z_17]);
702  }
703  }
704  }
705  }
706  if (omit_outside > 0)
707  printf("%d source space points omitted because they are outside the inner skull surface.\n",
708  omit_outside);
709  if (omit > 0)
710  printf("%d source space points omitted because of the %6.1f-mm distance limit.\n",
711  omit,1000*a->limit);
712  a->stat = OK;
713  return NULL;
714 }
715 
716 //=============================================================================================================
717 
718 int 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? */
719 /*
720  * Remove all source space points closer to the surface than a given limit
721  */
722 {
723  MneSurfaceOld* surf = NULL;
724  int k;
725  int nproc = QThread::idealThreadCount();
726  FilterThreadArg* a;
727 
728  if (!bemfile)
729  return OK;
730 
731  if ((surf = MneSurfaceOld::read_bem_surface(bemfile,FIFFV_BEM_SURF_ID_BRAIN,FALSE,NULL)) == NULL) {
732  qCritical("BEM model does not have the inner skull triangulation!");
733  return FAIL;
734  }
735  /*
736  * How close are the source points to the surface?
737  */
738  printf("Source spaces are in ");
739  if (spaces[0]->coord_frame == FIFFV_COORD_HEAD)
740  printf("head coordinates.\n");
741  else if (spaces[0]->coord_frame == FIFFV_COORD_MRI)
742  printf("MRI coordinates.\n");
743  else
744  printf("unknown (%d) coordinates.\n",spaces[0]->coord_frame);
745  printf("Checking that the sources are inside the inner skull ");
746  if (limit > 0.0)
747  printf("and at least %6.1f mm away",1000*limit);
748  printf(" (will take a few...)\n");
749  if (nproc < 2 || nspace == 1 || !use_threads) {
750  /*
751  * This is the conventional calculation
752  */
753  for (k = 0; k < nspace; k++) {
754  a = new FilterThreadArg();
755  a->s = spaces[k];
756  a->mri_head_t = mri_head_t;
757  a->surf = surf;
758  a->limit = limit;
759  a->filtered = filtered;
760  filter_source_space(a);
761  if(a)
762  delete a;
763  rearrange_source_space(spaces[k]);
764  }
765  }
766  else {
767  /*
768  * Calculate all (both) source spaces simultaneously
769  */
770  QList<FilterThreadArg*> args;//filterThreadArg *args = MALLOC_17(nspace,filterThreadArg);
771 
772  for (k = 0; k < nspace; k++) {
773  a = new FilterThreadArg();
774  a->s = spaces[k];
775  a->mri_head_t = mri_head_t;
776  a->surf = surf;
777  a->limit = limit;
778  a->filtered = filtered;
779  args.append(a);
780  }
781  /*
782  * Ready to start the threads & Wait for them to complete
783  */
784  QtConcurrent::blockingMap(args, filter_source_space);
785 
786  for (k = 0; k < nspace; k++) {
787  rearrange_source_space(spaces[k]);
788  if(args[k])
789  delete args[k];
790  }
791  }
792  if(surf)
793  delete surf;
794  printf("Thank you for waiting.\n\n");
795 
796  return OK;
797 }
798 
799 //=============================================================================================================
800 
801 MneSourceSpaceOld* MneSurfaceOrVolume::make_volume_source_space(MneSurfaceOld* surf, float grid, float exclude, float mindist)
802 /*
803  * Make a source space which covers the volume bounded by surf
804  */
805 {
806  float min[3],max[3],cm[3];
807  int minn[3],maxn[3];
808  float *node,maxdist,dist,diff[3];
809  int k,c;
810  MneSourceSpaceOld* sp = NULL;
811  int np,nplane,nrow;
812  int *neigh,nneigh;
813  int x,y,z;
814  /*
815  * Figure out the grid size
816  */
817  cm[X_17] = cm[Y_17] = cm[Z_17] = 0.0;
818  node = surf->rr[0];
819  for (c = 0; c < 3; c++)
820  min[c] = max[c] = node[c];
821 
822  for (k = 0; k < surf->np; k++) {
823  node = surf->rr[k];
824  for (c = 0; c < 3; c++) {
825  cm[c] += node[c];
826  if (node[c] < min[c])
827  min[c] = node[c];
828  if (node[c] > max[c])
829  max[c] = node[c];
830  }
831  }
832  for (c = 0; c < 3; c++)
833  cm[c] = cm[c]/surf->np;
834  /*
835  * Define the sphere which fits the surface
836  */
837  maxdist = 0.0;
838  for (k = 0; k < surf->np; k++) {
839  VEC_DIFF_17(cm,surf->rr[k],diff);
840  dist = VEC_LEN_17(diff);
841  if (dist > maxdist)
842  maxdist = dist;
843  }
844  printf("Surface CM = (%6.1f %6.1f %6.1f) mm\n",
845  1000*cm[X_17], 1000*cm[Y_17], 1000*cm[Z_17]);
846  printf("Surface fits inside a sphere with radius %6.1f mm\n",1000*maxdist);
847  printf("Surface extent:\n"
848  "\tx = %6.1f ... %6.1f mm\n"
849  "\ty = %6.1f ... %6.1f mm\n"
850  "\tz = %6.1f ... %6.1f mm\n",
851  1000*min[X_17],1000*max[X_17],
852  1000*min[Y_17],1000*max[Y_17],
853  1000*min[Z_17],1000*max[Z_17]);
854  for (c = 0; c < 3; c++) {
855  if (max[c] > 0)
856  maxn[c] = floor(std::fabs(max[c])/grid)+1;
857  else
858  maxn[c] = -floor(std::fabs(max[c])/grid)-1;
859  if (min[c] > 0)
860  minn[c] = floor(std::fabs(min[c])/grid)+1;
861  else
862  minn[c] = -floor(std::fabs(min[c])/grid)-1;
863  }
864  printf("Grid extent:\n"
865  "\tx = %6.1f ... %6.1f mm\n"
866  "\ty = %6.1f ... %6.1f mm\n"
867  "\tz = %6.1f ... %6.1f mm\n",
868  1000*(minn[X_17]*grid),1000*(maxn[X_17]*grid),
869  1000*(minn[Y_17]*grid),1000*(maxn[Y_17]*grid),
870  1000*(minn[Z_17]*grid),1000*(maxn[Z_17]*grid));
871  /*
872  * Now make the initial grid
873  */
874  np = 1;
875  for (c = 0; c < 3; c++)
876  np = np*(maxn[c]-minn[c]+1);
877  nplane = (maxn[X_17]-minn[X_17]+1)*(maxn[Y_17]-minn[Y_17]+1);
878  nrow = (maxn[X_17]-minn[X_17]+1);
879  sp = MneSurfaceOrVolume::mne_new_source_space(np);
880  sp->type = MNE_SOURCE_SPACE_VOLUME;
881  sp->nneighbor_vert = MALLOC_17(sp->np,int);
882  sp->neighbor_vert = MALLOC_17(sp->np,int *);
883  for (k = 0; k < sp->np; k++) {
884  sp->inuse[k] = TRUE;
885  sp->vertno[k] = k;
886  sp->nn[k][X_17] = sp->nn[k][Y_17] = 0.0; /* Source orientation is immaterial */
887  sp->nn[k][Z_17] = 1.0;
888  sp->neighbor_vert[k] = neigh = MALLOC_17(NNEIGHBORS,int);
889  sp->nneighbor_vert[k] = nneigh = NNEIGHBORS;
890  for (c = 0; c < nneigh; c++)
891  neigh[c] = -1;
892  sp->nuse++;
893  }
894  for (k = 0, z = minn[Z_17]; z <= maxn[Z_17]; z++) {
895  for (y = minn[Y_17]; y <= maxn[Y_17]; y++) {
896  for (x = minn[X_17]; x <= maxn[X_17]; x++, k++) {
897  sp->rr[k][X_17] = x*grid;
898  sp->rr[k][Y_17] = y*grid;
899  sp->rr[k][Z_17] = z*grid;
900  /*
901  * Figure out the neighborhood:
902  * 6-neighborhood first
903  */
904  neigh = sp->neighbor_vert[k];
905  if (z > minn[Z_17])
906  neigh[0] = k - nplane;
907  if (x < maxn[X_17])
908  neigh[1] = k + 1;
909  if (y < maxn[Y_17])
910  neigh[2] = k + nrow;
911  if (x > minn[X_17])
912  neigh[3] = k - 1;
913  if (y > minn[Y_17])
914  neigh[4] = k - nrow;
915  if (z < maxn[Z_17])
916  neigh[5] = k + nplane;
917  /*
918  * Then the rest to complete the 26-neighborhood
919  * First the plane below
920  */
921  if (z > minn[Z_17]) {
922  if (x < maxn[X_17]) {
923  neigh[6] = k + 1 - nplane;
924  if (y < maxn[Y_17])
925  neigh[7] = k + 1 + nrow - nplane;
926  }
927  if (y < maxn[Y_17])
928  neigh[8] = k + nrow - nplane;
929  if (x > minn[X_17]) {
930  if (y < maxn[Y_17])
931  neigh[9] = k - 1 + nrow - nplane;
932  neigh[10] = k - 1 - nplane;
933  if (y > minn[Y_17])
934  neigh[11] = k - 1 - nrow - nplane;
935  }
936  if (y > minn[Y_17]) {
937  neigh[12] = k - nrow - nplane;
938  if (x < maxn[X_17])
939  neigh[13] = k + 1 - nrow - nplane;
940  }
941  }
942  /*
943  * Then the same plane
944  */
945  if (x < maxn[X_17] && y < maxn[Y_17])
946  neigh[14] = k + 1 + nrow;
947  if (x > minn[X_17]) {
948  if (y < maxn[Y_17])
949  neigh[15] = k - 1 + nrow;
950  if (y > minn[Y_17])
951  neigh[16] = k - 1 - nrow;
952  }
953  if (y > minn[Y_17] && x < maxn[X_17])
954  neigh[17] = k + 1 - nrow - nplane;
955  /*
956  * Finally one plane above
957  */
958  if (z < maxn[Z_17]) {
959  if (x < maxn[X_17]) {
960  neigh[18] = k + 1 + nplane;
961  if (y < maxn[Y_17])
962  neigh[19] = k + 1 + nrow + nplane;
963  }
964  if (y < maxn[Y_17])
965  neigh[20] = k + nrow + nplane;
966  if (x > minn[X_17]) {
967  if (y < maxn[Y_17])
968  neigh[21] = k - 1 + nrow + nplane;
969  neigh[22] = k - 1 + nplane;
970  if (y > minn[Y_17])
971  neigh[23] = k - 1 - nrow + nplane;
972  }
973  if (y > minn[Y_17]) {
974  neigh[24] = k - nrow + nplane;
975  if (x < maxn[X_17])
976  neigh[25] = k + 1 - nrow + nplane;
977  }
978  }
979  }
980  }
981  }
982  printf("%d sources before omitting any.\n",sp->nuse);
983  /*
984  * Exclude infeasible points
985  */
986  for (k = 0; k < sp->np; k++) {
987  VEC_DIFF_17(cm,sp->rr[k],diff);
988  dist = VEC_LEN_17(diff);
989  if (dist < exclude || dist > maxdist) {
990  sp->inuse[k] = FALSE;
991  sp->nuse--;
992  }
993  }
994  printf("%d sources after omitting infeasible sources.\n",sp->nuse);
995  if (mne_filter_source_spaces(surf,mindist,NULL,&sp,1,NULL) != OK)
996  goto bad;
997  printf("%d sources remaining after excluding the sources outside the surface and less than %6.1f mm inside.\n",sp->nuse,1000*mindist);
998  /*
999  * Omit unused vertices from the neighborhoods
1000  */
1001  printf("Adjusting the neighborhood info...");
1002  for (k = 0; k < sp->np; k++) {
1003  neigh = sp->neighbor_vert[k];
1004  nneigh = sp->nneighbor_vert[k];
1005  if (sp->inuse[k]) {
1006  for (c = 0; c < nneigh; c++)
1007  if (!sp->inuse[neigh[c]])
1008  neigh[c] = -1;
1009  }
1010  else {
1011  for (c = 0; c < nneigh; c++)
1012  neigh[c] = -1;
1013  }
1014  }
1015  printf("[done]\n");
1016  /*
1017  * Set up the volume data (needed for creating the interpolation matrix)
1018  */
1019  {
1020  float r0[3],voxel_size[3],x_ras[3],y_ras[3],z_ras[3];
1021  int width,height,depth;
1022 
1023  r0[X_17] = minn[X_17]*grid;
1024  r0[Y_17] = minn[Y_17]*grid;
1025  r0[Z_17] = minn[Z_17]*grid;
1026 
1027  voxel_size[0] = grid;
1028  voxel_size[1] = grid;
1029  voxel_size[2] = grid;
1030 
1031  width = (maxn[X_17]-minn[X_17]+1);
1032  height = (maxn[Y_17]-minn[Y_17]+1);
1033  depth = (maxn[Z_17]-minn[Z_17]+1);
1034 
1035  for (k = 0; k < 3; k++)
1036  x_ras[k] = y_ras[k] = z_ras[k] = 0.0;
1037 
1038  x_ras[0] = 1.0;
1039  y_ras[1] = 1.0;
1040  z_ras[2] = 1.0;
1041 
1042  if ((sp->voxel_surf_RAS_t = make_voxel_ras_trans(r0,x_ras,y_ras,z_ras,voxel_size)) == NULL)
1043  goto bad;
1044 
1045  sp->vol_dims[0] = width;
1046  sp->vol_dims[1] = height;
1047  sp->vol_dims[2] = depth;
1048  VEC_COPY_17(sp->voxel_size,voxel_size);
1049  }
1050 
1051  return sp;
1052 
1053 bad : {
1054  if(sp)
1055  delete sp;
1056  return NULL;
1057  }
1058 }
1059 
1060 //=============================================================================================================
1061 
1062 MneSourceSpaceOld* MneSurfaceOrVolume::mne_new_source_space(int np)
1063 /*
1064  * Create a new source space and all associated data
1065  */
1066 {
1067  MneSourceSpaceOld* res = new MneSourceSpaceOld();
1068  res->np = np;
1069  if (np > 0) {
1070  res->rr = ALLOC_CMATRIX_17(np,3);
1071  res->nn = ALLOC_CMATRIX_17(np,3);
1072  res->inuse = ALLOC_INT_17(np);
1073  res->vertno = ALLOC_INT_17(np);
1074  }
1075  else {
1076  res->rr = NULL;
1077  res->nn = NULL;
1078  res->inuse = NULL;
1079  res->vertno = NULL;
1080  }
1081  res->nuse = 0;
1082  res->ntri = 0;
1083  res->tris = NULL;
1084  res->itris = NULL;
1085  res->tot_area = 0.0;
1086 
1087  res->nuse_tri = 0;
1088  res->use_tris = NULL;
1089  res->use_itris = NULL;
1090 
1091  res->neighbor_tri = NULL;
1092  res->nneighbor_tri = NULL;
1093  res->curv = NULL;
1094  res->val = NULL;
1095 
1096  res->neighbor_vert = NULL;
1097  res->nneighbor_vert = NULL;
1098  res->vert_dist = NULL;
1099 
1100  res->coord_frame = FIFFV_COORD_MRI;
1101  res->id = FIFFV_MNE_SURF_UNKNOWN;
1102  res->subject.clear();
1103  res->type = FIFFV_MNE_SPACE_SURFACE;
1104 
1105  res->nearest = NULL;
1106  res->patches = NULL;
1107  res->npatch = 0;
1108 
1109  res->dist = NULL;
1110  res->dist_limit = -1.0;
1111 
1112  res->voxel_surf_RAS_t = NULL;
1113  res->vol_dims[0] = res->vol_dims[1] = res->vol_dims[2] = 0;
1114 
1115  res->MRI_volume.clear();
1116  res->MRI_surf_RAS_RAS_t = NULL;
1117  res->MRI_voxel_surf_RAS_t = NULL;
1118  res->MRI_vol_dims[0] = res->MRI_vol_dims[1] = res->MRI_vol_dims[2] = 0;
1119  res->interpolator = NULL;
1120 
1121  res->vol_geom = NULL;
1122  res->mgh_tags = NULL;
1123  res->user_data = NULL;
1124  res->user_data_free = NULL;
1125 
1126  res->cm[X_17] = res->cm[Y_17] = res->cm[Z_17] = 0.0;
1127 
1128  return res;
1129 }
1130 
1131 //=============================================================================================================
1132 
1133 MneSurfaceOld* MneSurfaceOrVolume::read_bem_surface(const QString &name, int which, int add_geometry, float *sigmap) /* Conductivity? */
1134 {
1135  return read_bem_surface(name,which,add_geometry,sigmap,true);
1136 }
1137 
1138 //=============================================================================================================
1139 
1140 MneSurfaceOld* MneSurfaceOrVolume::mne_read_bem_surface2(char *name, int which, int add_geometry, float *sigmap)
1141 {
1142  return read_bem_surface(name,which,add_geometry,sigmap,FALSE);
1143 }
1144 
1145 //=============================================================================================================
1146 
1147 MneSurfaceOld* MneSurfaceOrVolume::read_bem_surface(const QString &name, int which, int add_geometry, float *sigmap, bool check_too_many_neighbors)
1148 /*
1149  * Read a Neuromag-style BEM surface description
1150  */
1151 {
1152  QFile file(name);
1153  FiffStream::SPtr stream(new FiffStream(&file));
1154 
1155  QList<FiffDirNode::SPtr> surfs;
1156  QList<FiffDirNode::SPtr> bems;
1157  FiffDirNode::SPtr node;
1158  FiffTag::SPtr t_pTag;
1159 
1160  int id = -1;
1161  float **nodes = NULL;
1162  float **node_normals = NULL;
1163  int **triangles = NULL;
1164  int nnode,ntri;
1165  MneSurfaceOld* s = NULL;
1166  int k;
1167  int coord_frame = FIFFV_COORD_MRI;
1168  float sigma = -1.0;
1169  MatrixXf tmp_nodes;
1170  MatrixXi tmp_triangles;
1171 
1172  if(!stream->open())
1173  goto bad;
1174  /*
1175  * Check for the existence of BEM coord frame
1176  */
1177  bems = stream->dirtree()->dir_tree_find(FIFFB_BEM);
1178  if (bems.size() > 0) {
1179  node = bems[0];
1180  if (node->find_tag(stream, FIFF_BEM_COORD_FRAME, t_pTag)) {
1181  coord_frame = *t_pTag->toInt();
1182  }
1183  }
1184  surfs = stream->dirtree()->dir_tree_find(FIFFB_BEM_SURF);
1185  if (surfs.size() == 0) {
1186  printf ("No BEM surfaces found in %s",name.toUtf8().constData());
1187  goto bad;
1188  }
1189  if (which >= 0) {
1190  for (k = 0; k < surfs.size(); ++k) {
1191  node = surfs[k];
1192  /*
1193  * Read the data from this node
1194  */
1195  if (node->find_tag(stream, FIFF_BEM_SURF_ID, t_pTag)) {
1196  id = *t_pTag->toInt();
1197  if (id == which)
1198  break;
1199  }
1200  }
1201  if (id != which) {
1202  printf("Desired surface not found in %s",name.toUtf8().constData());
1203  goto bad;
1204  }
1205  }
1206  else
1207  node = surfs[0];
1208  /*
1209  * Get the compulsory tags
1210  */
1211  if (!node->find_tag(stream, FIFF_BEM_SURF_NNODE, t_pTag))
1212  goto bad;
1213  nnode = *t_pTag->toInt();
1214 
1215  if (!node->find_tag(stream, FIFF_BEM_SURF_NTRI, t_pTag))
1216  goto bad;
1217  ntri = *t_pTag->toInt();
1218 
1219  if (!node->find_tag(stream, FIFF_BEM_SURF_NODES, t_pTag))
1220  goto bad;
1221  tmp_nodes = t_pTag->toFloatMatrix().transpose();
1222  nodes = ALLOC_CMATRIX_17(tmp_nodes.rows(),tmp_nodes.cols());
1223  fromFloatEigenMatrix_17(tmp_nodes, nodes);
1224 
1225  if (node->find_tag(stream, FIFF_BEM_SURF_NORMALS, t_pTag)) {\
1226  MatrixXf tmp_node_normals = t_pTag->toFloatMatrix().transpose();
1227  node_normals = ALLOC_CMATRIX_17(tmp_node_normals.rows(),tmp_node_normals.cols());
1228  fromFloatEigenMatrix_17(tmp_node_normals, node_normals);
1229  }
1230 
1231  if (!node->find_tag(stream, FIFF_BEM_SURF_TRIANGLES, t_pTag))
1232  goto bad;
1233  tmp_triangles = t_pTag->toIntMatrix().transpose();
1234  triangles = (int **)malloc(tmp_triangles.rows() * sizeof(int *));
1235  for (int i = 0; i < tmp_triangles.rows(); ++i)
1236  triangles[i] = (int *)malloc(tmp_triangles.cols() * sizeof(int));
1237  fromIntEigenMatrix_17(tmp_triangles, triangles);
1238 
1239  if (node->find_tag(stream, FIFF_MNE_COORD_FRAME, t_pTag)) {
1240  coord_frame = *t_pTag->toInt();
1241  }
1242  else if (node->find_tag(stream, FIFF_BEM_COORD_FRAME, t_pTag)) {
1243  coord_frame = *t_pTag->toInt();
1244  }
1245  if (node->find_tag(stream, FIFF_BEM_SIGMA, t_pTag)) {
1246  sigma = *t_pTag->toFloat();
1247  }
1248 
1249  stream->close();
1250 
1251  s = (MneSurfaceOld*)mne_new_source_space(0);
1252  for (k = 0; k < ntri; k++) {
1253  triangles[k][0]--;
1254  triangles[k][1]--;
1255  triangles[k][2]--;
1256  }
1257  s->itris = triangles;
1258  s->id = which;
1259  s->coord_frame = coord_frame;
1260  s->rr = nodes; nodes = NULL;
1261  s->nn = node_normals; node_normals = NULL;
1262  s->ntri = ntri;
1263  s->np = nnode;
1264  s->curv = NULL;
1265  s->val = NULL;
1266 
1267  if (add_geometry) {
1268  if (check_too_many_neighbors) {
1269  if (mne_source_space_add_geometry_info((MneSourceSpaceOld*)s,!s->nn) != OK)
1270  goto bad;
1271  }
1272  else {
1273  if (mne_source_space_add_geometry_info2((MneSourceSpaceOld*)s,!s->nn) != OK)
1274  goto bad;
1275  }
1276  }
1277  else if (s->nn == NULL) { /* Normals only */
1278  if (mne_add_vertex_normals((MneSourceSpaceOld*)s) != OK)
1279  goto bad;
1280  }
1281  else
1282  mne_add_triangle_data((MneSourceSpaceOld*)s);
1283 
1284  s->nuse = s->np;
1285  s->inuse = MALLOC_17(s->np,int);
1286  s->vertno = MALLOC_17(s->np,int);
1287  for (k = 0; k < s->np; k++) {
1288  s->inuse[k] = TRUE;
1289  s->vertno[k] = k;
1290  }
1291  if (sigmap)
1292  *sigmap = sigma;
1293 
1294  return s;
1295 
1296 bad : {
1297  FREE_CMATRIX_17(nodes);
1298  FREE_CMATRIX_17(node_normals);
1299  FREE_ICMATRIX_17(triangles);
1300  stream->close();
1301  return NULL;
1302  }
1303 }
1304 
1305 //=============================================================================================================
1306 
1307 void MneSurfaceOrVolume::mne_triangle_coords(float *r, MneSurfaceOld* s, int tri, float *x, float *y, float *z)
1308 /*
1309  * Compute the coordinates of a point within a triangle
1310  */
1311 {
1312  double rr[3]; /* Vector from triangle corner #1 to r */
1313  double a,b,c,v1,v2,det;
1314  MneTriangle* this_tri;
1315 
1316  this_tri = s->tris+tri;
1317 
1318  VEC_DIFF_17(this_tri->r1,r,rr);
1319  *z = VEC_DOT_17(rr,this_tri->nn);
1320 
1321  a = VEC_DOT_17(this_tri->r12,this_tri->r12);
1322  b = VEC_DOT_17(this_tri->r13,this_tri->r13);
1323  c = VEC_DOT_17(this_tri->r12,this_tri->r13);
1324 
1325  v1 = VEC_DOT_17(rr,this_tri->r12);
1326  v2 = VEC_DOT_17(rr,this_tri->r13);
1327 
1328  det = a*b - c*c;
1329 
1330  *x = (b*v1 - c*v2)/det;
1331  *y = (a*v2 - c*v1)/det;
1332 
1333  return;
1334 }
1335 
1336 //=============================================================================================================
1337 
1338 int MneSurfaceOrVolume::nearest_triangle_point(float *r, MneSurfaceOld* s, void *user, int tri, float *x, float *y, float *z)
1339 /*
1340  * Find the nearest point from a triangle
1341  */
1342 {
1343 
1344  double p,q,p0,q0,t0;
1345  double rr[3]; /* Vector from triangle corner #1 to r */
1346  double a,b,c,v1,v2,det;
1347  double best,dist,dist0;
1348  MneProjData* pd = (MneProjData*)user;
1349  MneTriangle* this_tri;
1350 
1351  this_tri = s->tris+tri;
1352  VEC_DIFF_17(this_tri->r1,r,rr);
1353  dist = VEC_DOT_17(rr,this_tri->nn);
1354 
1355  if (pd) {
1356  if (!pd->act[tri])
1357  return FALSE;
1358  a = pd->a[tri];
1359  b = pd->b[tri];
1360  c = pd->c[tri];
1361  }
1362  else {
1363  a = VEC_DOT_17(this_tri->r12,this_tri->r12);
1364  b = VEC_DOT_17(this_tri->r13,this_tri->r13);
1365  c = VEC_DOT_17(this_tri->r12,this_tri->r13);
1366  }
1367 
1368  v1 = VEC_DOT_17(rr,this_tri->r12);
1369  v2 = VEC_DOT_17(rr,this_tri->r13);
1370 
1371  det = a*b - c*c;
1372 
1373  p = (b*v1 - c*v2)/det;
1374  q = (a*v2 - c*v1)/det;
1375  /*
1376  * If the point projects into the triangle we are done
1377  */
1378  if (p >= 0.0 && p <= 1.0 &&
1379  q >= 0.0 && q <= 1.0 &&
1380  q <= 1.0 - p) {
1381  *x = p;
1382  *y = q;
1383  *z = dist;
1384  return TRUE;
1385  }
1386  /*
1387  * Tough: must investigate the sides
1388  * We might do something intelligent here. However, for now it is ok
1389  * to do it in the hard way
1390  */
1391  /*
1392  * Side 1 -> 2
1393  */
1394  p0 = p + 0.5*(q * c)/a;
1395  if (p0 < 0.0)
1396  p0 = 0.0;
1397  else if (p0 > 1.0)
1398  p0 = 1.0;
1399  q0 = 0.0;
1400  dist0 = sqrt((p-p0)*(p-p0)*a +
1401  (q-q0)*(q-q0)*b +
1402  (p-p0)*(q-q0)*c +
1403  dist*dist);
1404  best = dist0;
1405  *x = p0;
1406  *y = q0;
1407  *z = dist0;
1408  /*
1409  * Side 2 -> 3
1410  */
1411  t0 = 0.5*((2.0*a-c)*(1.0-p) + (2.0*b-c)*q)/(a+b-c);
1412  if (t0 < 0.0)
1413  t0 = 0.0;
1414  else if (t0 > 1.0)
1415  t0 = 1.0;
1416  p0 = 1.0 - t0;
1417  q0 = t0;
1418  dist0 = sqrt((p-p0)*(p-p0)*a +
1419  (q-q0)*(q-q0)*b +
1420  (p-p0)*(q-q0)*c +
1421  dist*dist);
1422  if (dist0 < best) {
1423  best = dist0;
1424  *x = p0;
1425  *y = q0;
1426  *z = dist0;
1427  }
1428  /*
1429  * Side 1 -> 3
1430  */
1431  p0 = 0.0;
1432  q0 = q + 0.5*(p * c)/b;
1433  if (q0 < 0.0)
1434  q0 = 0.0;
1435  else if (q0 > 1.0)
1436  q0 = 1.0;
1437  dist0 = sqrt((p-p0)*(p-p0)*a +
1438  (q-q0)*(q-q0)*b +
1439  (p-p0)*(q-q0)*c +
1440  dist*dist);
1441  if (dist0 < best) {
1442  best = dist0;
1443  *x = p0;
1444  *y = q0;
1445  *z = dist0;
1446  }
1447  return TRUE;
1448 }
1449 
1450 //=============================================================================================================
1451 
1452 void MneSurfaceOrVolume::project_to_triangle(MneSurfaceOld* s, int tri, float p, float q, float *r)
1453 {
1454  int k;
1455  MneTriangle* this_tri;
1456 
1457  this_tri = s->tris+tri;
1458 
1459  for (k = 0; k < 3; k++)
1460  r[k] = this_tri->r1[k] + p*this_tri->r12[k] + q*this_tri->r13[k];
1461 
1462  return;
1463 }
1464 
1465 //=============================================================================================================
1466 
1467 int MneSurfaceOrVolume::mne_nearest_triangle_point(float *r, MneSurfaceOld* s, int tri, float *x, float *y, float *z)
1468 /*
1469  * This is for external use
1470  */
1471 {
1472  return nearest_triangle_point(r,s,NULL,tri,x,y,z);
1473 }
1474 
1475 //=============================================================================================================
1476 
1477 int MneSurfaceOrVolume::mne_project_to_surface(MneSurfaceOld* s, void *proj_data, float *r, int project_it, float *distp)
1478 /*
1479  * Project the point onto the closest point on the surface
1480  */
1481 {
1482  float dist; /* Distance to the triangle */
1483  float p,q; /* Coordinates on the triangle */
1484  float p0,q0,dist0;
1485  int best;
1486  int k;
1487 
1488  p0 = q0 = 0.0;
1489  dist0 = 0.0;
1490  for (best = -1, k = 0; k < s->ntri; k++) {
1491  if (nearest_triangle_point(r,s,proj_data,k,&p,&q,&dist)) {
1492  if (best < 0 || std::fabs(dist) < std::fabs(dist0)) {
1493  dist0 = dist;
1494  best = k;
1495  p0 = p;
1496  q0 = q;
1497  }
1498  }
1499  }
1500  if (best >= 0 && project_it)
1501  project_to_triangle(s,best,p0,q0,r);
1502  if (distp)
1503  *distp = dist0;
1504  return best;
1505 }
1506 
1507 //=============================================================================================================
1508 
1509 void MneSurfaceOrVolume::mne_project_to_triangle(MneSurfaceOld* s,
1510  int best,
1511  float *r,
1512  float *proj)
1513 /*
1514  * Project to a triangle provided that we know the best match already
1515  */
1516 {
1517  float p,q,dist;
1518 
1519  mne_nearest_triangle_point(r,s,best,&p,&q,&dist);
1520  project_to_triangle(s,best,p,q,proj);
1521 
1522  return;
1523 }
1524 
1525 //=============================================================================================================
1526 
1527 void MneSurfaceOrVolume::mne_find_closest_on_surface_approx(MneSurfaceOld* s, float **r, int np, int *nearest, float *dist, int nstep)
1528 /*
1529  * Find the closest triangle on the surface for each point and the distance to it
1530  * This uses the values in nearest as approximations of the closest triangle
1531  */
1532 {
1533  MneProjData* p = new MneProjData(s);
1534  int k,was;
1535  float mydist;
1536 
1537  printf("%s for %d points %d steps...",nearest[0] < 0 ? "Closest" : "Approx closest",np,nstep);
1538 
1539  for (k = 0; k < np; k++) {
1540  was = nearest[k];
1541  decide_search_restriction(s,p,nearest[k],nstep,r[k]);
1542  nearest[k] = mne_project_to_surface(s,p,r[k],0,dist ? dist+k : &mydist);
1543  if (nearest[k] < 0) {
1544  decide_search_restriction(s,p,-1,nstep,r[k]);
1545  nearest[k] = mne_project_to_surface(s,p,r[k],0,dist ? dist+k : &mydist);
1546  }
1547  }
1548 
1549  printf("[done]\n");
1550  delete p;
1551  return;
1552 }
1553 
1554 //=============================================================================================================
1555 
1556 void MneSurfaceOrVolume::decide_search_restriction(MneSurfaceOld* s,
1557  MneProjData* p,
1558  int approx_best, /* We know the best triangle approximately
1559  * already */
1560  int nstep,
1561  float *r)
1562 /*
1563  * Restrict the search only to feasible triangles
1564  */
1565 {
1566  int k;
1567  float diff[3],dist,mindist;
1568  int minvert;
1569 
1570  for (k = 0; k < s->ntri; k++)
1571  p->act[k] = FALSE;
1572 
1573  if (approx_best < 0) {
1574  /*
1575  * Search for the closest vertex
1576  */
1577  mindist = 1000.0;
1578  minvert = 0;
1579  for (k = 0; k < s->np; k++) {
1580  VEC_DIFF_17(r,s->rr[k],diff);
1581  dist = VEC_LEN_17(diff);
1582  if (dist < mindist && s->nneighbor_tri[k] > 0) {
1583  mindist = dist;
1584  minvert = k;
1585  }
1586  }
1587  }
1588  else {
1589  /*
1590  * Just use this triangle
1591  */
1592  MneTriangle* this_tri = NULL;
1593 
1594  this_tri = s->tris+approx_best;
1595  VEC_DIFF_17(r,this_tri->r1,diff);
1596  mindist = VEC_LEN_17(diff);
1597  minvert = this_tri->vert[0];
1598 
1599  VEC_DIFF_17(r,this_tri->r2,diff);
1600  dist = VEC_LEN_17(diff);
1601  if (dist < mindist) {
1602  mindist = dist;
1603  minvert = this_tri->vert[1];
1604  }
1605  VEC_DIFF_17(r,this_tri->r3,diff);
1606  dist = VEC_LEN_17(diff);
1607  if (dist < mindist) {
1608  mindist = dist;
1609  minvert = this_tri->vert[2];
1610  }
1611  }
1612  /*
1613  * Activate triangles in the neighborhood
1614  */
1615  activate_neighbors(s,minvert,p->act,nstep);
1616 
1617  for (k = 0, p->nactive = 0; k < s->ntri; k++)
1618  if (p->act[k])
1619  p->nactive++;
1620  return;
1621 }
1622 
1623 //=============================================================================================================
1624 
1625 void MneSurfaceOrVolume::activate_neighbors(MneSurfaceOld* s, int start, int *act, int nstep)
1626 /*
1627  * Blessed recursion...
1628  */
1629 {
1630  int k;
1631 
1632  if (nstep == 0)
1633  return;
1634 
1635  for (k = 0; k < s->nneighbor_tri[start]; k++)
1636  act[s->neighbor_tri[start][k]] = TRUE;
1637  for (k = 0; k < s->nneighbor_vert[start]; k++)
1638  activate_neighbors(s,s->neighbor_vert[start][k],act,nstep-1);
1639 
1640  return;
1641 }
1642 
1643 //=============================================================================================================
1644 
1645 int MneSurfaceOrVolume::mne_read_source_spaces(const QString &name, MneSourceSpaceOld* **spacesp, int *nspacep)
1646 /*
1647  * Read source spaces from a FIFF file
1648  */
1649 {
1650  QFile file(name);
1651  FiffStream::SPtr stream(new FiffStream(&file));
1652 
1653  int nspace = 0;
1654  MneSourceSpaceOld* *spaces = NULL;
1655  MneSourceSpaceOld* new_space = NULL;
1656  QList<FiffDirNode::SPtr> sources;
1657  FiffDirNode::SPtr node;
1658  FiffTag::SPtr t_pTag;
1659  int j,k,p,q;
1660  int ntri;
1661  int *nearest = NULL;
1662  float *nearest_dist = NULL;
1663  int *nneighbors = NULL;
1664  int *neighbors = NULL;
1665  int *vol_dims = NULL;
1666 
1667  if(!stream->open())
1668  goto bad;
1669 
1670  sources = stream->dirtree()->dir_tree_find(FIFFB_MNE_SOURCE_SPACE);
1671  if (sources.size() == 0) {
1672  printf("No source spaces available here");
1673  goto bad;
1674  }
1675  for (j = 0; j < sources.size(); j++) {
1676  new_space = MneSurfaceOrVolume::mne_new_source_space(0);
1677  node = sources[j];
1678  /*
1679  * Get the mandatory data first
1680  */
1681  if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NPOINTS, t_pTag)) {
1682  goto bad;
1683  }
1684  new_space->np = *t_pTag->toInt();
1685  if (new_space->np == 0) {
1686  printf("No points in this source space");
1687  goto bad;
1688  }
1689  if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_POINTS, t_pTag)) {
1690  goto bad;
1691  }
1692  MatrixXf tmp_rr = t_pTag->toFloatMatrix().transpose();
1693  new_space->rr = ALLOC_CMATRIX_17(tmp_rr.rows(),tmp_rr.cols());
1694  fromFloatEigenMatrix_17(tmp_rr,new_space->rr);
1695  if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NORMALS, t_pTag)) {
1696  goto bad;
1697  }
1698  MatrixXf tmp_nn = t_pTag->toFloatMatrix().transpose();
1699  new_space->nn = ALLOC_CMATRIX_17(tmp_nn.rows(),tmp_nn.cols());
1700  fromFloatEigenMatrix_17(tmp_nn,new_space->nn);
1701  if (!node->find_tag(stream, FIFF_MNE_COORD_FRAME, t_pTag)) {
1702  new_space->coord_frame = FIFFV_COORD_MRI;
1703  }
1704  else {
1705  new_space->coord_frame = *t_pTag->toInt();
1706  }
1707  if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_ID, t_pTag)) {
1708  new_space->id = *t_pTag->toInt();
1709  }
1710  if (node->find_tag(stream, FIFF_SUBJ_HIS_ID, t_pTag)) {
1711  new_space->subject = (char *)t_pTag->data();
1712  }
1713  if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_TYPE, t_pTag)) {
1714  new_space->type = *t_pTag->toInt();
1715  }
1716  ntri = 0;
1717  if (node->find_tag(stream, FIFF_BEM_SURF_NTRI, t_pTag)) {
1718  ntri = *t_pTag->toInt();
1719  }
1720  else if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NTRI, t_pTag)) {
1721  ntri = *t_pTag->toInt();
1722  }
1723  if (ntri > 0) {
1724  int **itris = NULL;
1725 
1726  if (!node->find_tag(stream, FIFF_BEM_SURF_TRIANGLES, t_pTag)) {
1727  if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_TRIANGLES, t_pTag)) {
1728  goto bad;
1729  }
1730  }
1731 
1732  MatrixXi tmp_itris = t_pTag->toIntMatrix().transpose();
1733  itris = (int **)malloc(tmp_itris.rows() * sizeof(int *));
1734  for (int i = 0; i < tmp_itris.rows(); ++i)
1735  itris[i] = (int *)malloc(tmp_itris.cols() * sizeof(int));
1736  fromIntEigenMatrix_17(tmp_itris, itris);
1737 
1738  for (p = 0; p < ntri; p++) { /* Adjust the numbering */
1739  itris[p][X_17]--;
1740  itris[p][Y_17]--;
1741  itris[p][Z_17]--;
1742  }
1743  new_space->itris = itris; itris = NULL;
1744  new_space->ntri = ntri;
1745  }
1746  if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NUSE, t_pTag)) {
1747  if (new_space->type == FIFFV_MNE_SPACE_VOLUME) {
1748  /*
1749  * Use all
1750  */
1751  new_space->nuse = new_space->np;
1752  new_space->inuse = MALLOC_17(new_space->nuse,int);
1753  new_space->vertno = MALLOC_17(new_space->nuse,int);
1754  for (k = 0; k < new_space->nuse; k++) {
1755  new_space->inuse[k] = TRUE;
1756  new_space->vertno[k] = k;
1757  }
1758  }
1759  else {
1760  /*
1761  * None in use
1762  * NOTE: The consequences of this change have to be evaluated carefully
1763  */
1764  new_space->nuse = 0;
1765  new_space->inuse = MALLOC_17(new_space->np,int);
1766  new_space->vertno = NULL;
1767  for (k = 0; k < new_space->np; k++)
1768  new_space->inuse[k] = FALSE;
1769  }
1770  }
1771  else {
1772  new_space->nuse = *t_pTag->toInt();
1773  if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_SELECTION, t_pTag)) {
1774  goto bad;
1775  }
1776 
1777  qDebug() << "ToDo: Check whether new_space->inuse contains the right stuff!!! - use VectorXi instead";
1778  // new_space->inuse = t_pTag->toInt();
1779  new_space->inuse = MALLOC_17(new_space->np,int); //DEBUG
1780  if (new_space->nuse > 0) {
1781  new_space->vertno = MALLOC_17(new_space->nuse,int);
1782  for (k = 0, p = 0; k < new_space->np; k++) {
1783  new_space->inuse[k] = t_pTag->toInt()[k]; //DEBUG
1784  if (new_space->inuse[k])
1785  new_space->vertno[p++] = k;
1786  }
1787  }
1788  else {
1789  FREE_17(new_space->vertno);
1790  new_space->vertno = NULL;
1791  }
1792  /*
1793  * Selection triangulation
1794  */
1795  ntri = 0;
1796  if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NUSE_TRI, t_pTag)) {
1797  ntri = *t_pTag->toInt();
1798  }
1799  if (ntri > 0) {
1800  int **itris = NULL;
1801 
1802  if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_USE_TRIANGLES, t_pTag)) {
1803  goto bad;
1804  }
1805 
1806  MatrixXi tmp_itris = t_pTag->toIntMatrix().transpose();
1807  itris = (int **)malloc(tmp_itris.rows() * sizeof(int *));
1808  for (int i = 0; i < tmp_itris.rows(); ++i)
1809  itris[i] = (int *)malloc(tmp_itris.cols() * sizeof(int));
1810  fromIntEigenMatrix_17(tmp_itris, itris);
1811  for (p = 0; p < ntri; p++) { /* Adjust the numbering */
1812  itris[p][X_17]--;
1813  itris[p][Y_17]--;
1814  itris[p][Z_17]--;
1815  }
1816  new_space->use_itris = itris; itris = NULL;
1817  new_space->nuse_tri = ntri;
1818  }
1819  /*
1820  * The patch information becomes relevant here
1821  */
1822  if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NEAREST, t_pTag)) {
1823  nearest = t_pTag->toInt();
1824  new_space->nearest = MALLOC_17(new_space->np,MneNearest);
1825  for (k = 0; k < new_space->np; k++) {
1826  new_space->nearest[k].vert = k;
1827  new_space->nearest[k].nearest = nearest[k];
1828  new_space->nearest[k].patch = NULL;
1829  }
1830 
1831  if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NEAREST_DIST, t_pTag)) {
1832  goto bad;
1833  }
1834  qDebug() << "ToDo: Check whether nearest_dist contains the right stuff!!! - use VectorXf instead";
1835  nearest_dist = t_pTag->toFloat();
1836  for (k = 0; k < new_space->np; k++) {
1837  new_space->nearest[k].dist = nearest_dist[k];
1838  }
1839  // FREE_17(nearest); nearest = NULL;
1840  // FREE_17(nearest_dist); nearest_dist = NULL;
1841  }
1842  /*
1843  * We may have the distance matrix
1844  */
1845  if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_DIST_LIMIT, t_pTag)) {
1846  new_space->dist_limit = *t_pTag->toFloat();
1847  if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_DIST, t_pTag)) {
1848  // SparseMatrix<double> tmpSparse = t_pTag->toSparseFloatMatrix();
1849  FiffSparseMatrix* dist = FiffSparseMatrix::fiff_get_float_sparse_matrix(t_pTag);
1850  new_space->dist = dist->mne_add_upper_triangle_rcs();
1851  delete dist;
1852  if (!new_space->dist) {
1853  goto bad;
1854  }
1855  }
1856  else
1857  new_space->dist_limit = 0.0;
1858  }
1859  }
1860  /*
1861  * For volume source spaces we might have the neighborhood information
1862  */
1863  if (new_space->type == FIFFV_MNE_SPACE_VOLUME) {
1864  int ntot,nvert,ntot_count,nneigh;
1865  int *neigh;
1866 
1867  nneighbors = neighbors = NULL;
1868  ntot = nvert = 0;
1869  if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NEIGHBORS, t_pTag)) {
1870  qDebug() << "ToDo: Check whether neighbors contains the right stuff!!! - use VectorXi instead";
1871  neighbors = t_pTag->toInt();
1872  ntot = t_pTag->size()/sizeof(fiff_int_t);
1873  }
1874  if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NNEIGHBORS, t_pTag)) {
1875  qDebug() << "ToDo: Check whether nneighbors contains the right stuff!!! - use VectorXi instead";
1876  nneighbors = t_pTag->toInt();
1877  nvert = t_pTag->size()/sizeof(fiff_int_t);
1878  }
1879  if (neighbors && nneighbors) {
1880  if (nvert != new_space->np) {
1881  printf("Inconsistent neighborhood data in file.");
1882  goto bad;
1883  }
1884  for (k = 0, ntot_count = 0; k < nvert; k++)
1885  ntot_count += nneighbors[k];
1886  if (ntot_count != ntot) {
1887  printf("Inconsistent neighborhood data in file.");
1888  goto bad;
1889  }
1890  new_space->nneighbor_vert = MALLOC_17(nvert,int);
1891  new_space->neighbor_vert = MALLOC_17(nvert,int *);
1892  for (k = 0, q = 0; k < nvert; k++) {
1893  new_space->nneighbor_vert[k] = nneigh = nneighbors[k];
1894  new_space->neighbor_vert[k] = neigh = MALLOC_17(nneigh,int);
1895  for (p = 0; p < nneigh; p++,q++)
1896  neigh[p] = neighbors[q];
1897  }
1898  }
1899  FREE_17(neighbors);
1900  FREE_17(nneighbors);
1901  nneighbors = neighbors = NULL;
1902  /*
1903  * There might be a coordinate transformation and dimensions
1904  */
1905  new_space->voxel_surf_RAS_t = FiffCoordTransOld::mne_read_transform_from_node(stream, node, FIFFV_MNE_COORD_MRI_VOXEL, FIFFV_MNE_COORD_SURFACE_RAS);
1906  if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_VOXEL_DIMS, t_pTag)) {
1907  qDebug() << "ToDo: Check whether vol_dims contains the right stuff!!! - use VectorXi instead";
1908  vol_dims = t_pTag->toInt();
1909  }
1910  if (vol_dims)
1911  VEC_COPY_17(new_space->vol_dims,vol_dims);
1912  {
1913  QList<FiffDirNode::SPtr> mris = node->dir_tree_find(FIFFB_MNE_PARENT_MRI_FILE);
1914 
1915  if (mris.size() == 0) { /* The old way */
1916  new_space->MRI_surf_RAS_RAS_t = FiffCoordTransOld::mne_read_transform_from_node(stream, node, FIFFV_MNE_COORD_SURFACE_RAS, FIFFV_MNE_COORD_RAS);
1917  if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_MRI_FILE, t_pTag)) {
1918  qDebug() << "ToDo: Check whether new_space->MRI_volume contains the right stuff!!! - use QString instead";
1919  new_space->MRI_volume = (char *)t_pTag->data();
1920  }
1921  if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_INTERPOLATOR, t_pTag)) {
1922  new_space->interpolator = FiffSparseMatrix::fiff_get_float_sparse_matrix(t_pTag);
1923  }
1924  }
1925  else {
1926  if (node->find_tag(stream, FIFF_MNE_FILE_NAME, t_pTag)) {
1927  new_space->MRI_volume = (char *)t_pTag->data();
1928  }
1929  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);
1930  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);
1931 
1932  if (mris[0]->find_tag(stream, FIFF_MNE_SOURCE_SPACE_INTERPOLATOR, t_pTag)) {
1933  new_space->interpolator = FiffSparseMatrix::fiff_get_float_sparse_matrix(t_pTag);
1934  }
1935  if (mris[0]->find_tag(stream, FIFF_MRI_WIDTH, t_pTag)) {
1936  new_space->MRI_vol_dims[0] = *t_pTag->toInt();
1937  }
1938  if (mris[0]->find_tag(stream, FIFF_MRI_HEIGHT, t_pTag)) {
1939  new_space->MRI_vol_dims[1] = *t_pTag->toInt();
1940  }
1941  if (mris[0]->find_tag(stream, FIFF_MRI_DEPTH, t_pTag)) {
1942  new_space->MRI_vol_dims[2] = *t_pTag->toInt();
1943  }
1944  }
1945  }
1946  }
1947  mne_add_triangle_data(new_space);
1948  spaces = REALLOC_17(spaces,nspace+1,MneSourceSpaceOld*);
1949  spaces[nspace++] = new_space;
1950  new_space = NULL;
1951  }
1952  stream->close();
1953 
1954  *spacesp = spaces;
1955  *nspacep = nspace;
1956 
1957  return FIFF_OK;
1958 
1959 bad : {
1960  stream->close();
1961  delete new_space;
1962  for (k = 0; k < nspace; k++)
1963  delete spaces[k];
1964  FREE_17(spaces);
1965  FREE_17(nearest);
1966  FREE_17(nearest_dist);
1967  FREE_17(neighbors);
1968  FREE_17(nneighbors);
1969  FREE_17(vol_dims);
1970 
1971  return FIFF_FAIL;
1972  }
1973 }
1974 
1975 //=============================================================================================================
1976 
1977 void MneSurfaceOrVolume::mne_source_space_update_inuse(MneSourceSpaceOld* s, int *new_inuse)
1978 /*
1979  * Update the active vertices
1980  */
1981 {
1982  int k,p,nuse;
1983 
1984  if (!s)
1985  return;
1986 
1987  FREE_17(s->inuse); s->inuse = new_inuse;
1988 
1989  for (k = 0, nuse = 0; k < s->np; k++)
1990  if (s->inuse[k])
1991  nuse++;
1992 
1993  s->nuse = nuse;
1994  if (s->nuse > 0) {
1995  s->vertno = REALLOC_17(s->vertno,s->nuse,int);
1996  for (k = 0, p = 0; k < s->np; k++)
1997  if (s->inuse[k])
1998  s->vertno[p++] = k;
1999  }
2000  else {
2001  FREE_17(s->vertno);
2002  s->vertno = NULL;
2003  }
2004  return;
2005 }
2006 
2007 //=============================================================================================================
2008 
2009 int MneSurfaceOrVolume::mne_is_left_hemi_source_space(MneSourceSpaceOld* s)
2010 /*
2011  * Left or right hemisphere?
2012  */
2013 {
2014  int k;
2015  float xave;
2016 
2017  for (k = 0, xave = 0.0; k < s->np; k++)
2018  xave += s->rr[k][0];
2019  if (xave < 0.0)
2020  return TRUE;
2021  else
2022  return FALSE;
2023 }
2024 
2025 //=============================================================================================================
2026 
2027 int MneSurfaceOrVolume::mne_transform_source_space(MneSourceSpaceOld* ss, FiffCoordTransOld *t)
2028 /*
2029  * Transform source space data into another coordinate frame
2030  */
2031 {
2032  int k;
2033  if (ss == NULL)
2034  return OK;
2035  if (ss->coord_frame == t->to)
2036  return OK;
2037  if (ss->coord_frame != t->from) {
2038  printf("Coordinate transformation does not match with the source space coordinate system.");
2039  return FAIL;
2040  }
2041  for (k = 0; k < ss->np; k++) {
2042  FiffCoordTransOld::fiff_coord_trans(ss->rr[k],t,FIFFV_MOVE);
2043  FiffCoordTransOld::fiff_coord_trans(ss->nn[k],t,FIFFV_NO_MOVE);
2044  }
2045  if (ss->tris) {
2046  for (k = 0; k < ss->ntri; k++)
2047  FiffCoordTransOld::fiff_coord_trans(ss->tris[k].nn,t,FIFFV_NO_MOVE);
2048  }
2049  ss->coord_frame = t->to;
2050  return OK;
2051 }
2052 
2053 //=============================================================================================================
2054 
2055 int MneSurfaceOrVolume::mne_transform_source_spaces_to(int coord_frame, FiffCoordTransOld *t, MneSourceSpaceOld* *spaces, int nspace)
2056 /*
2057  * Facilitate the transformation of the source spaces
2058  */
2059 {
2060  MneSourceSpaceOld* s;
2061  int k;
2062  FiffCoordTransOld* my_t;
2063 
2064  for (k = 0; k < nspace; k++) {
2065  s = spaces[k];
2066  if (s->coord_frame != coord_frame) {
2067  if (t) {
2068  if (s->coord_frame == t->from && t->to == coord_frame) {
2069  if (mne_transform_source_space(s,t) != OK)
2070  return FAIL;
2071  }
2072  else if (s->coord_frame == t->to && t->from == coord_frame) {
2073  my_t = t->fiff_invert_transform();
2074  if (mne_transform_source_space(s,my_t) != OK) {
2075  FREE_17(my_t);
2076  return FAIL;
2077  }
2078  FREE_17(my_t);
2079  }
2080  else {
2081  printf("Could not transform a source space because of transformation incompatibility.");
2082  return FAIL;
2083  }
2084  }
2085  else {
2086  printf("Could not transform a source space because of missing coordinate transformation.");
2087  return FAIL;
2088  }
2089  }
2090  }
2091  return OK;
2092 }
2093 
2094 //=============================================================================================================
2095 
2096 void MneSurfaceOrVolume::enable_all_sources(MneSourceSpaceOld* s)
2097 {
2098  int k;
2099  for (k = 0; k < s->np; k++)
2100  s->inuse[k] = TRUE;
2101  s->nuse = s->np;
2102  return;
2103 }
2104 
2105 //=============================================================================================================
2106 
2107 #define LH_LABEL_TAG "-lh.label"
2108 #define RH_LABEL_TAG "-rh.label"
2109 
2110 int MneSurfaceOrVolume::restrict_sources_to_labels(MneSourceSpaceOld* *spaces, int nspace, const QStringList& labels, int nlabel)
2111 /*
2112  * Pick only sources within a label
2113  */
2114 {
2115  MneSourceSpaceOld* lh = NULL;
2116  MneSourceSpaceOld* rh = NULL;
2117  MneSourceSpaceOld* sp;
2118  int *lh_inuse = NULL;
2119  int *rh_inuse = NULL;
2120  int *sel = NULL;
2121  int nsel;
2122  int *inuse;
2123  int k,p;
2124 
2125  if (nlabel == 0)
2126  return OK;
2127 
2128  for (k = 0; k < nspace; k++) {
2129  if (mne_is_left_hemi_source_space(spaces[k])) {
2130  lh = spaces[k];
2131  FREE_17(lh_inuse);
2132  lh_inuse = MALLOC_17(lh->np,int);
2133  for (p = 0; p < lh->np; p++)
2134  lh_inuse[p] = 0;
2135  }
2136  else {
2137  rh = spaces[k];
2138  FREE_17(rh_inuse);
2139  rh_inuse = MALLOC_17(rh->np,int);
2140  for (p = 0; p < rh->np; p++)
2141  rh_inuse[p] = 0;
2142  }
2143  }
2144  /*
2145  * Go through each label file
2146  */
2147  for (k = 0; k < nlabel; k++) {
2148  /*
2149  * Which hemi?
2150  */
2151  if (labels[k].contains(LH_LABEL_TAG)){ //strstr(labels[k],LH_LABEL_TAG) != NULL) {
2152  sp = lh;
2153  inuse = lh_inuse;
2154  }
2155  else if (labels[k].contains(RH_LABEL_TAG)){ //strstr(labels[k],RH_LABEL_TAG) != NULL) {
2156  sp = rh;
2157  inuse = rh_inuse;
2158  }
2159  else {
2160  printf("\tWarning: cannot assign label file %s to a hemisphere.\n",labels[k].toUtf8().constData());
2161  continue;
2162  }
2163  if (sp) {
2164  if (mne_read_label(labels[k],NULL,&sel,&nsel) == FAIL)
2165  goto bad;
2166  for (p = 0; p < nsel; p++) {
2167  if (sel[p] >= 0 && sel[p] < sp->np)
2168  inuse[sel[p]] = sp->inuse[sel[p]];
2169  else
2170  printf("vertex number out of range in %s (%d vs %d)\n",
2171  labels[k].toUtf8().constData(),sel[p],sp->np);
2172  }
2173  FREE_17(sel); sel = NULL;
2174  printf("Processed label file %s\n",labels[k].toUtf8().constData());
2175  }
2176  }
2177  mne_source_space_update_inuse(lh,lh_inuse);
2178  mne_source_space_update_inuse(rh,rh_inuse);
2179  return OK;
2180 
2181 bad : {
2182  FREE_17(lh_inuse);
2183  FREE_17(rh_inuse);
2184  FREE_17(sel);
2185  return FAIL;
2186  }
2187 }
2188 
2189 //=============================================================================================================
2190 
2191 int MneSurfaceOrVolume::mne_find_sources_in_label(char *label, MneSourceSpaceOld* s, int off, int **selp, int *nselp) /* How many selected? */
2192 /*
2193  * Find the source points within a label
2194  */
2195 {
2196  FILE *in = NULL;
2197  int res = FAIL;
2198 
2199  int nsel = 0;
2200  int *sel = NULL;
2201 
2202  int k,p,pp,nlabel,q;
2203  char c;
2204  float fdum;
2205  /*
2206  * Read the label file
2207  */
2208  if ((in = fopen(label,"r")) == NULL) {
2209  qCritical(label);//err_set_sys_error(label);
2210  goto out;
2211  }
2212  c = fgetc(in);
2213  if (c !='#') {
2214  qCritical("Label file does not start correctly.");
2215  goto out;
2216  }
2217  for (c = fgetc(in); c != '\n'; c = fgetc(in))
2218  ;
2219  if (fscanf(in,"%d",&nlabel) != 1) {
2220  qCritical("Could not read the number of labelled points.");
2221  goto out;
2222  }
2223 #ifdef DEBUG
2224  printf("\t%d points in label %s\n",nlabel,label);
2225 #endif
2226  for (k = 0; k < nlabel; k++) {
2227  if (fscanf(in,"%d %g %g %g %g",&p,
2228  &fdum, &fdum, &fdum, &fdum) != 5) {
2229  qCritical("Could not read label point # %d",k+1);
2230  goto out;
2231  }
2232  if (p < 0 || p >= s->np) {
2233  qCritical("Source index out of range %d (range 0..%d)\n",p,s->np-1);
2234  goto out;
2235  }
2236  if (s->inuse[p]) {
2237  for (pp = 0, q = 0; pp < p; pp++) {
2238  if (s->inuse[pp])
2239  q++;
2240  }
2241  sel = REALLOC_17(sel,nsel+1,int);
2242  sel[nsel++] = q + off;
2243  }
2244  }
2245  *nselp = nsel;
2246  *selp = sel;
2247  res = OK;
2248 
2249 out : {
2250  if (in)
2251  fclose(in);
2252  if (res != OK) {
2253  FREE_17(sel);
2254  *selp = NULL;
2255  *nselp = 0;
2256  }
2257  return res;
2258  }
2259 }
2260 
2261 //=============================================================================================================
2262 
2263 int MneSurfaceOrVolume::mne_read_label(const QString& label, char **commentp, int **selp, int *nselp) /* How many? */
2264 /*
2265  * Find the source points within a label
2266  */
2267 {
2268  FILE *in = NULL;
2269  int res = FAIL;
2270 
2271  int nsel = 0;
2272  int *sel = NULL;
2273 
2274  int k,p,nlabel;
2275  char c;
2276  char *comment = NULL;
2277  float fdum;
2278  /*
2279  * Read the label file
2280  */
2281  if ((in = fopen(label.toUtf8().constData(),"r")) == NULL) {
2282  qCritical() << label;//err_set_sys_error(label);
2283  goto out;
2284  }
2285  for (k = 0; k < 2; k++) {
2286  rewind(in);
2287  c = fgetc(in);
2288  if (c !='#') {
2289  qCritical("Label file does not start correctly.");
2290  goto out;
2291  }
2292  for (c = fgetc(in); c == ' ' && c != '\n'; c = fgetc(in))
2293  ;
2294  ungetc(c,in);
2295  if (k == 0) {
2296  for (c = fgetc(in), p = 0; c != '\n'; c = fgetc(in), p++)
2297  ;
2298  }
2299  else {
2300  for (c = fgetc(in); c != '\n'; c = fgetc(in))
2301  *comment++ = c;
2302  *comment = '\0';
2303  }
2304  if (!commentp)
2305  break;
2306  if (p == 0) {
2307  *commentp = NULL;
2308  break;
2309  }
2310  if (k == 0)
2311  *commentp = comment = MALLOC_17(p+1,char);
2312  }
2313  if (fscanf(in,"%d",&nlabel) != 1) {
2314  qCritical("Could not read the number of labelled points.");
2315  goto out;
2316  }
2317  for (k = 0; k < nlabel; k++) {
2318  if (fscanf(in,"%d %g %g %g %g",&p,
2319  &fdum, &fdum, &fdum, &fdum) != 5) {
2320  qCritical("Could not read label point # %d",k+1);
2321  goto out;
2322  }
2323  sel = REALLOC_17(sel,nsel+1,int);
2324  sel[nsel++] = p;
2325  }
2326  *nselp = nsel;
2327  *selp = sel;
2328  res = OK;
2329 
2330 out : {
2331  if (in)
2332  fclose(in);
2333  if (res != OK) {
2334  FREE_17(sel);
2335  *selp = NULL;
2336  *nselp = 0;
2337  }
2338  return res;
2339  }
2340 }
2341 
2342 //=============================================================================================================
2343 
2344 int MneSurfaceOrVolume::mne_write_label(char *label, char *comment, int *sel, int nsel, float **rr) /* Locations of the nodes in MRI coords */
2345 /*
2346  * Find the source points within a label
2347  */
2348 {
2349  FILE *out = NULL;
2350  int res = FAIL;
2351  int k;
2352  float fdum = 0.0;
2353  /*
2354  * Read the label file
2355  */
2356  if ((out = fopen(label,"w")) == NULL) {
2357  qCritical(label);//err_set_sys_error(label);
2358  goto out;
2359  }
2360  if (comment == NULL)
2361  fprintf(out,"# Label file created by the MNE software\n");
2362  else
2363  fprintf(out,"# %s\n",comment);
2364 
2365  fprintf(out,"%d\n",nsel);
2366  if (rr != NULL)
2367  for (k = 0; k < nsel; k++)
2368  fprintf(out,"%d %.2f %.2f %.2f %g\n",sel[k],
2369  1000*rr[sel[k]][0],
2370  1000*rr[sel[k]][1],
2371  1000*rr[sel[k]][2],fdum);
2372  else
2373  for (k = 0; k < nsel; k++)
2374  fprintf(out,"%d %.2f %.2f %.2f %g\n",sel[k],fdum,fdum,fdum,fdum);
2375  res = OK;
2376 
2377 out : {
2378  if (out)
2379  fclose(out);
2380  if (res != OK)
2381  unlink(label);
2382  return res;
2383  }
2384 }
2385 
2386 //=============================================================================================================
2387 
2388 void MneSurfaceOrVolume::mne_add_triangle_data(MneSourceSpaceOld* s)
2389 /*
2390  * Add the triangle data structures
2391  */
2392 {
2393  int k;
2394  MneTriangle* tri;
2395 
2396  if (!s || s->type != MNE_SOURCE_SPACE_SURFACE)
2397  return;
2398 
2399  FREE_17(s->tris); s->tris = NULL;
2400  FREE_17(s->use_tris); s->use_tris = NULL;
2401  /*
2402  * Add information for the complete triangulation
2403  */
2404  if (s->itris && s->ntri > 0) {
2405  s->tris = MALLOC_17(s->ntri,MneTriangle);
2406  s->tot_area = 0.0;
2407  for (k = 0, tri = s->tris; k < s->ntri; k++, tri++) {
2408  tri->vert = s->itris[k];
2409  tri->r1 = s->rr[tri->vert[0]];
2410  tri->r2 = s->rr[tri->vert[1]];
2411  tri->r3 = s->rr[tri->vert[2]];
2412  MneTriangle::add_triangle_data(tri);
2413  s->tot_area += tri->area;
2414  }
2415 #ifdef TRIANGLE_SIZE_WARNING
2416  for (k = 0, tri = s->tris; k < s->ntri; k++, tri++)
2417  if (tri->area < 1e-5*s->tot_area/s->ntri)
2418  printf("Warning: Triangle area is only %g um^2 (%.5f %% of expected average)\n",
2419  1e12*tri->area,100*s->ntri*tri->area/s->tot_area);
2420 #endif
2421  }
2422 #ifdef DEBUG
2423  printf("\ttotal area = %-.1f cm^2\n",1e4*s->tot_area);
2424 #endif
2425  /*
2426  * Add information for the selected subset if applicable
2427  */
2428  if (s->use_itris && s->nuse_tri > 0) {
2429  s->use_tris = MALLOC_17(s->nuse_tri,MneTriangle);
2430  for (k = 0, tri = s->use_tris; k < s->nuse_tri; k++, tri++) {
2431  tri->vert = s->use_itris[k];
2432  tri->r1 = s->rr[tri->vert[0]];
2433  tri->r2 = s->rr[tri->vert[1]];
2434  tri->r3 = s->rr[tri->vert[2]];
2435  MneTriangle::add_triangle_data(tri);
2436  }
2437  }
2438  return;
2439 }
2440 
2441 //=============================================================================================================
2442 
2443 void MneSurfaceOrVolume::mne_compute_cm(float **rr, int np, float *cm)
2444 /*
2445  * Compute the center of mass of a set of points
2446  */
2447 {
2448  int q;
2449  cm[0] = cm[1] = cm[2] = 0.0;
2450  for (q = 0; q < np; q++) {
2451  cm[0] += rr[q][0];
2452  cm[1] += rr[q][1];
2453  cm[2] += rr[q][2];
2454  }
2455  if (np > 0) {
2456  cm[0] = cm[0]/np;
2457  cm[1] = cm[1]/np;
2458  cm[2] = cm[2]/np;
2459  }
2460  return;
2461 }
2462 
2463 //=============================================================================================================
2464 
2465 void MneSurfaceOrVolume::mne_compute_surface_cm(MneSurfaceOld *s)
2466 /*
2467  * Compute the center of mass of a surface
2468  */
2469 {
2470  if (!s)
2471  return;
2472 
2473  mne_compute_cm(s->rr,s->np,s->cm);
2474  return;
2475 }
2476 
2477 //=============================================================================================================
2478 
2479 void MneSurfaceOrVolume::calculate_vertex_distances(MneSourceSpaceOld* s)
2480 {
2481  int k,p,ndist;
2482  float *dist,diff[3];
2483  int *neigh, nneigh;
2484 
2485  if (!s->neighbor_vert || !s->nneighbor_vert)
2486  return;
2487 
2488  if (s->vert_dist) {
2489  for (k = 0; k < s->np; k++)
2490  FREE_17(s->vert_dist[k]);
2491  FREE_17(s->vert_dist);
2492  }
2493  s->vert_dist = MALLOC_17(s->np,float *);
2494  printf("\tDistances between neighboring vertices...");
2495  for (k = 0, ndist = 0; k < s->np; k++) {
2496  s->vert_dist[k] = dist = MALLOC_17(s->nneighbor_vert[k],float);
2497  neigh = s->neighbor_vert[k];
2498  nneigh = s->nneighbor_vert[k];
2499  for (p = 0; p < nneigh; p++) {
2500  if (neigh[p] >= 0) {
2501  VEC_DIFF_17(s->rr[k],s->rr[neigh[p]],diff);
2502  dist[p] = VEC_LEN_17(diff);
2503  }
2504  else
2505  dist[p] = -1.0;
2506  ndist++;
2507  }
2508  }
2509  printf("[%d distances done]\n",ndist);
2510  return;
2511 }
2512 
2513 //=============================================================================================================
2514 
2515 int MneSurfaceOrVolume::mne_add_vertex_normals(MneSourceSpaceOld* s)
2516 {
2517  int k,c,p;
2518  int *ii;
2519  float w,size;
2520  MneTriangle* tri;
2521 
2522  if (!s || s->type != MNE_SOURCE_SPACE_SURFACE)
2523  return OK;
2524  /*
2525  * Reallocate the stuff and initialize
2526  */
2527  FREE_CMATRIX_17(s->nn);
2528  s->nn = ALLOC_CMATRIX_17(s->np,3);
2529 
2530  for (k = 0; k < s->np; k++) {
2531  s->nn[k][X_17] = s->nn[k][Y_17] = s->nn[k][Z_17] = 0.0;
2532  }
2533  /*
2534  * One pass through the triangles will do it
2535  */
2536  MneSurfaceOrVolume::mne_add_triangle_data(s);
2537  for (p = 0, tri = s->tris; p < s->ntri; p++, tri++) {
2538  ii = tri->vert;
2539  w = 1.0; /* This should be related to the triangle size */
2540  /*
2541  * Then the vertex normals
2542  */
2543  for (k = 0; k < 3; k++)
2544  for (c = 0; c < 3; c++)
2545  s->nn[ii[k]][c] += w*tri->nn[c];
2546  }
2547  for (k = 0; k < s->np; k++) {
2548  size = VEC_LEN_17(s->nn[k]);
2549  if (size > 0.0)
2550  for (c = 0; c < 3; c++)
2551  s->nn[k][c] = s->nn[k][c]/size;
2552  }
2553  mne_compute_surface_cm((MneSurfaceOld*)s);
2554  return OK;
2555 }
2556 
2557 //=============================================================================================================
2558 
2559 int MneSurfaceOrVolume::add_geometry_info(MneSourceSpaceOld* s, int do_normals, int *border, int check_too_many_neighbors)
2560 /*
2561  * Add vertex normals and neighbourhood information
2562  */
2563 {
2564  int k,c,p,q;
2565  int vert;
2566  int *ii;
2567  int *neighbors,nneighbors;
2568  float w,size;
2569  int found;
2570  int nfix_distinct,nfix_no_neighbors,nfix_defect;
2571  MneTriangle* tri;
2572 
2573  if (!s)
2574  return OK;
2575 
2576  if (s->type == MNE_SOURCE_SPACE_VOLUME) {
2577  calculate_vertex_distances(s);
2578  return OK;
2579  }
2580  if (s->type != MNE_SOURCE_SPACE_SURFACE)
2581  return OK;
2582  /*
2583  * Reallocate the stuff and initialize
2584  */
2585  if (do_normals) {
2586  FREE_CMATRIX_17(s->nn);
2587  s->nn = ALLOC_CMATRIX_17(s->np,3);
2588  }
2589  if (s->neighbor_tri) {
2590  for (k = 0; k < s->np; k++)
2591  FREE_17(s->neighbor_tri[k]);
2592  FREE_17(s->neighbor_tri);
2593  }
2594  FREE_17(s->nneighbor_tri);
2595  s->neighbor_tri = MALLOC_17(s->np,int *);
2596  s->nneighbor_tri = MALLOC_17(s->np,int);
2597 
2598  for (k = 0; k < s->np; k++) {
2599  s->neighbor_tri[k] = NULL;
2600  s->nneighbor_tri[k] = 0;
2601  if (do_normals)
2602  s->nn[k][X_17] = s->nn[k][Y_17] = s->nn[k][Z_17] = 0.0;
2603  }
2604  /*
2605  * One pass through the triangles will do it
2606  */
2607  mne_add_triangle_data(s);
2608  for (p = 0, tri = s->tris; p < s->ntri; p++, tri++)
2609  if (tri->area == 0)
2610  printf("\tWarning : zero size triangle # %d\n",p);
2611  printf("\tTriangle ");
2612  if (do_normals)
2613  printf("and vertex ");
2614  printf("normals and neighboring triangles...");
2615  for (p = 0, tri = s->tris; p < s->ntri; p++, tri++) {
2616  ii = tri->vert;
2617  w = 1.0; /* This should be related to the triangle size */
2618  for (k = 0; k < 3; k++) {
2619  /*
2620  * Then the vertex normals
2621  */
2622  if (do_normals)
2623  for (c = 0; c < 3; c++)
2624  s->nn[ii[k]][c] += w*tri->nn[c];
2625  /*
2626  * Add to the list of neighbors
2627  */
2628  s->neighbor_tri[ii[k]] = REALLOC_17(s->neighbor_tri[ii[k]],
2629  s->nneighbor_tri[ii[k]]+1,int);
2630  s->neighbor_tri[ii[k]][s->nneighbor_tri[ii[k]]] = p;
2631  s->nneighbor_tri[ii[k]]++;
2632  }
2633  }
2634  nfix_no_neighbors = 0;
2635  nfix_defect = 0;
2636  for (k = 0; k < s->np; k++) {
2637  if (s->nneighbor_tri[k] <= 0) {
2638  if (!border || !border[k]) {
2639 #ifdef STRICT_ERROR
2640  err_printf_set_error("Vertex %d does not have any neighboring triangles!",k);
2641  return FAIL;
2642 #else
2643 #ifdef REPORT_WARNINGS
2644  printf("Warning: Vertex %d does not have any neighboring triangles!\n",k);
2645 #endif
2646 #endif
2647  nfix_no_neighbors++;
2648  }
2649  }
2650  else if (s->nneighbor_tri[k] < 3 && !border) {
2651 #ifdef REPORT_WARNINGS
2652  printf("\n\tTopological defect: Vertex %d has only %d neighboring triangle%s Vertex omitted.\n\t",
2653  k,s->nneighbor_tri[k],s->nneighbor_tri[k] > 1 ? "s." : ".");
2654 #endif
2655  nfix_defect++;
2656  s->nneighbor_tri[k] = 0;
2657  FREE_17(s->neighbor_tri[k]);
2658  s->neighbor_tri[k] = NULL;
2659  }
2660  }
2661  /*
2662  * Scale the vertex normals to unit length
2663  */
2664  for (k = 0; k < s->np; k++)
2665  if (s->nneighbor_tri[k] > 0) {
2666  size = VEC_LEN_17(s->nn[k]);
2667  if (size > 0.0)
2668  for (c = 0; c < 3; c++)
2669  s->nn[k][c] = s->nn[k][c]/size;
2670  }
2671  printf("[done]\n");
2672  /*
2673  * Determine the neighboring vertices
2674  */
2675  printf("\tVertex neighbors...");
2676  if (s->neighbor_vert) {
2677  for (k = 0; k < s->np; k++)
2678  FREE_17(s->neighbor_vert[k]);
2679  FREE_17(s->neighbor_vert);
2680  }
2681  FREE_17(s->nneighbor_vert);
2682  s->neighbor_vert = MALLOC_17(s->np,int *);
2683  s->nneighbor_vert = MALLOC_17(s->np,int);
2684  /*
2685  * We know the number of neighbors beforehand
2686  */
2687  if (border) {
2688  for (k = 0; k < s->np; k++) {
2689  if (s->nneighbor_tri[k] > 0) {
2690  if (border[k]) {
2691  s->neighbor_vert[k] = MALLOC_17(s->nneighbor_tri[k]+1,int);
2692  s->nneighbor_vert[k] = s->nneighbor_tri[k]+1;
2693  }
2694  else {
2695  s->neighbor_vert[k] = MALLOC_17(s->nneighbor_tri[k],int);
2696  s->nneighbor_vert[k] = s->nneighbor_tri[k];
2697  }
2698  }
2699  else {
2700  s->neighbor_vert[k] = NULL;
2701  s->nneighbor_vert[k] = 0;
2702  }
2703  }
2704  }
2705  else {
2706  for (k = 0; k < s->np; k++) {
2707  if (s->nneighbor_tri[k] > 0) {
2708  s->neighbor_vert[k] = MALLOC_17(s->nneighbor_tri[k],int);
2709  s->nneighbor_vert[k] = s->nneighbor_tri[k];
2710  }
2711  else {
2712  s->neighbor_vert[k] = NULL;
2713  s->nneighbor_vert[k] = 0;
2714  }
2715  }
2716  }
2717  nfix_distinct = 0;
2718  for (k = 0; k < s->np; k++) {
2719  neighbors = s->neighbor_vert[k];
2720  nneighbors = 0;
2721  for (p = 0; p < s->nneighbor_tri[k]; p++) {
2722  /*
2723  * Fit in the other vertices of the neighboring triangle
2724  */
2725  for (c = 0; c < 3; c++) {
2726  vert = s->tris[s->neighbor_tri[k][p]].vert[c];
2727  if (vert != k) {
2728  for (q = 0, found = FALSE; q < nneighbors; q++) {
2729  if (neighbors[q] == vert) {
2730  found = TRUE;
2731  break;
2732  }
2733  }
2734  if (!found) {
2735  if (nneighbors < s->nneighbor_vert[k])
2736  neighbors[nneighbors++] = vert;
2737  else if (!border || !border[k]) {
2738  if (check_too_many_neighbors) {
2739  printf("Too many neighbors for vertex %d.",k);
2740  return FAIL;
2741  }
2742  else
2743  printf("\tWarning: Too many neighbors for vertex %d\n",k);
2744  }
2745  }
2746  }
2747  }
2748  }
2749  if (nneighbors != s->nneighbor_vert[k]) {
2750 #ifdef REPORT_WARNINGS
2751  printf("\n\tIncorrect number of distinct neighbors for vertex %d (%d instead of %d) [fixed].",
2752  k,nneighbors,s->nneighbor_vert[k]);
2753 #endif
2754  nfix_distinct++;
2755  s->nneighbor_vert[k] = nneighbors;
2756  }
2757  }
2758  printf("[done]\n");
2759  /*
2760  * Distance calculation follows
2761  */
2762  calculate_vertex_distances(s);
2763  mne_compute_surface_cm((MneSurfaceOld*)s);
2764  /*
2765  * Summarize the defects
2766  */
2767  if (nfix_defect > 0)
2768  printf("\tWarning: %d topological defects were fixed.\n",nfix_defect);
2769  if (nfix_distinct > 0)
2770  printf("\tWarning: %d vertices had incorrect number of distinct neighbors (fixed).\n",nfix_distinct);
2771  if (nfix_no_neighbors > 0)
2772  printf("\tWarning: %d vertices did not have any neighboring triangles (fixed)\n",nfix_no_neighbors);
2773 #ifdef DEBUG
2774  for (k = 0; k < s->np; k++) {
2775  if (s->nneighbor_vert[k] <= 0)
2776  printf("No neighbors for vertex %d\n",k);
2777  if (s->nneighbor_tri[k] <= 0)
2778  printf("No neighbor tris for vertex %d\n",k);
2779  }
2780 #endif
2781  return OK;
2782 }
2783 
2784 //=============================================================================================================
2785 
2786 int MneSurfaceOrVolume::mne_source_space_add_geometry_info(MneSourceSpaceOld* s, int do_normals)
2787 {
2788  return add_geometry_info(s,do_normals,NULL,TRUE);
2789 }
2790 
2791 //=============================================================================================================
2792 
2793 int MneSurfaceOrVolume::mne_source_space_add_geometry_info2(MneSourceSpaceOld* s, int do_normals)
2794 
2795 {
2796  return add_geometry_info(s,do_normals,NULL,FALSE);
2797 }
2798 
2799 //=============================================================================================================
2800 
2801 int MneSurfaceOrVolume::mne_label_area(char *label, MneSourceSpaceOld* s, float *areap) /* Return the area here */
2802 /*
2803  * Calculate the area of the label
2804  */
2805 {
2806  int *sel = NULL;
2807  int nsel = 0;
2808  float area;
2809  int k,q,nneigh,*neigh;
2810 
2811  if (!s) {
2812  qCritical("Source space not specified for mne_label_area");
2813  goto bad;
2814  }
2815  if (mne_read_label(label,NULL,&sel,&nsel))
2816  goto bad;
2817 
2818  area = 0.0;
2819  for (k = 0; k < nsel; k++) {
2820  if (sel[k] < 0 || sel[k] >= s->np) {
2821  qCritical("Label vertex index out of range in mne_label_area");
2822  goto bad;
2823  }
2824  nneigh = s->nneighbor_tri[sel[k]];
2825  neigh = s->neighbor_tri[sel[k]];
2826  for (q = 0; q < nneigh; q++)
2827  area += s->tris[neigh[q]].area/3.0;
2828  }
2829  FREE_17(sel);
2830  *areap = area;
2831  return OK;
2832 
2833 bad : {
2834  FREE_17(sel);
2835  return FAIL;
2836  }
2837 }
2838 
2839 //=============================================================================================================
2840 
2841 // Align the MEG fiducials to the MRI fiducials
2842 int MneSurfaceOrVolume::align_fiducials(FiffDigitizerData* head_dig,
2843  FiffDigitizerData* mri_dig,
2844  MneMshDisplaySurface* head_surf,
2845  int niter,
2846  int scale_head,
2847  float omit_dist,
2848  float *scales)
2849 
2850 {
2851  float *head_fid[3],*mri_fid[3],**fid;
2852  int j,k;
2853  FiffDigPoint p;
2854  FiffDigitizerData* dig = NULL;
2855  float nasion_weight = 5.0;
2856 
2857 
2858 
2859  if (!head_dig) {
2860  qCritical("MEG head coordinate system digitizer data not available");
2861  goto bad;
2862  }
2863  if (!mri_dig) {
2864  qCritical("MRI coordinate system digitizer data not available");
2865  goto bad;
2866  }
2867 
2868  for (j = 0; j < 2; j++) {
2869  dig = j == 0 ? head_dig : mri_dig;
2870  fid = j == 0 ? head_fid : mri_fid;
2871 
2872  for (k = 0; k < 3; k++) {
2873  fid[k] = NULL;
2874  }
2875 
2876  for (k = 0; k < dig->npoint; k++) {
2877  p = dig->points[k];
2878  if (p.kind == FIFFV_POINT_CARDINAL) {
2879  if (p.ident == FIFFV_POINT_LPA) {
2880  fid[0] = dig->points[k].r;
2881  }
2882  else if (p.ident == FIFFV_POINT_NASION) {
2883  fid[1] = dig->points[k].r;
2884  }
2885  else if (p.ident == FIFFV_POINT_RPA) {
2886  fid[2] = dig->points[k].r;
2887  }
2888  }
2889  }
2890  }
2891 
2892  for (k = 0; k < 3; k++) {
2893  if (!head_fid[k]) {
2894  qCritical("Some of the MEG fiducials were missing");
2895  goto bad;
2896  }
2897 
2898  if (!mri_fid[k]) {
2899  qCritical("Some of the MRI fiducials were missing");
2900  goto bad;
2901  }
2902  }
2903 
2904  if (scale_head) {
2905  get_head_scale(head_dig,mri_fid,head_surf,scales);
2906  printf("xscale = %.3f yscale = %.3f zscale = %.3f\n",scales[0],scales[1],scales[2]);
2907 
2908  for (j = 0; j < 3; j++)
2909  for (k = 0; k < 3; k++)
2910  mri_fid[j][k] = mri_fid[j][k]*scales[k];
2911 
2912  scale_display_surface(head_surf,scales);
2913  }
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 
2942 bad :
2943  return FAIL;
2944 }
2945 
2946 //=============================================================================================================
2947 
2948 // Simple head size fit
2949 void 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) == FAIL){
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) == FAIL) {
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 
3006 out : {
3007  FREE_17(dig_rr);
3008  FREE_17(head_rr);
3009  return;
3010  }
3011 }
3012 
3013 //=============================================================================================================
3014 
3015 int 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 
3047 void 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 
3136 int 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 
3213 out : {
3214  FREE_CMATRIX_17(rr_head);
3215  FREE_CMATRIX_17(rr_mri);
3216  FREE_17(w);
3217  return res;
3218  }
3219 }
3220 
3221 //=============================================================================================================
3222 
3223 float 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 
3242 void 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 
3265 void 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 
3280 char * 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 
3314 MneSourceSpaceOld* 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 
3322 MneSourceSpaceOld* 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 
3393 bad : {
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 
3405 int 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(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 
3613 bad : {
3614  if (fp)
3615  fclose(fp);
3616  FREE_CMATRIX_17(vert);
3617  FREE_ICMATRIX_17(tri);
3618  return FAIL;
3619  }
3620 }
3621 
3622 //=============================================================================================================
3623 
3624 int 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(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 
3710 bad : {
3711  if (fp)
3712  fclose(fp);
3713  FREE_17(curvs);
3714  return FAIL;
3715  }
3716 }
3717 
3718 //=============================================================================================================
3719 
3720 int 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 
3742 int 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 
3754 MneVolGeom* 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 
3775 MneVolGeom* 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 
3788 int 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 
3811 int 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 
3852 int 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 
3898 MneMghTagGroup* 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 
3915 MneVolGeom* 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  return vg;
3999 }
4000 
4001 //=============================================================================================================
4002 
4003 int MneSurfaceOrVolume::mne_read_int3(FILE *in, int *ival)
4004 /*
4005  * Read the strange 3-byte integer
4006  */
4007 {
4008  unsigned int s = 0;
4009 
4010  if (fread (&s,3,1,in) != 1) {
4011  if (ferror(in))
4012  qCritical("mne_read_int3");
4013  else
4014  qCritical("mne_read_int3 could not read data");
4015  return FAIL;
4016  }
4017  s = (unsigned int)UTILSLIB::IOUtils::swap_int(s);
4018  *ival = ((s >> 8) & 0xffffff);
4019  return OK;
4020 }
4021 
4022 //=============================================================================================================
4023 
4024 int MneSurfaceOrVolume::mne_read_int(FILE *in, qint32 *ival)
4025 /*
4026  * Read a 32-bit integer
4027  */
4028 {
4029  qint32 s ;
4030  if (fread (&s,sizeof(qint32),1,in) != 1) {
4031  if (ferror(in))
4032  qCritical("mne_read_int");
4033  else
4034  qCritical("mne_read_int could not read data");
4035  return FAIL;
4036  }
4037  *ival = UTILSLIB::IOUtils::swap_int(s);
4038  return OK;
4039 }
4040 
4041 //=============================================================================================================
4042 
4043 int MneSurfaceOrVolume::mne_read_int2(FILE *in, int *ival)
4044 /*
4045  * Read int from short
4046  */
4047 {
4048  short s ;
4049  if (fread (&s,sizeof(short),1,in) != 1) {
4050  if (ferror(in))
4051  qCritical("mne_read_int2");
4052  else
4053  qCritical("mne_read_int2 could not read data");
4054  return FAIL;
4055  }
4056  *ival = UTILSLIB::IOUtils::swap_short(s);
4057  return OK;
4058 }
4059 
4060 //=============================================================================================================
4061 
4062 int MneSurfaceOrVolume::mne_read_float(FILE *in, float *fval)
4063 /*
4064  * Read float
4065  */
4066 {
4067  float f ;
4068  if (fread (&f,sizeof(float),1,in) != 1) {
4069  if (ferror(in))
4070  qCritical("mne_read_float");
4071  else
4072  qCritical("mne_read_float could not read data");
4073  return FAIL;
4074  }
4075  *fval = UTILSLIB::IOUtils::swap_float(f);
4076  return OK;
4077 }
4078 
4079 //=============================================================================================================
4080 
4081 int MneSurfaceOrVolume::mne_read_long(FILE *in, long long *lval)
4082 /*
4083  * Read a 64-bit integer
4084  */
4085 {
4086  long long s ;
4087  if (fread (&s,sizeof(long long),1,in) != 1) {
4088  if (ferror(in))
4089  qCritical("mne_read_long");
4090  else
4091  qCritical("mne_read_long could not read data");
4092  return FAIL;
4093  }
4094  *lval = UTILSLIB::IOUtils::swap_long(s);
4095  return OK;
4096 }
4097 
4098 //=============================================================================================================
4099 
4100 char *MneSurfaceOrVolume::mne_strdup(const char *s)
4101 {
4102  char *res;
4103  if (s == NULL)
4104  return NULL;
4105  res = MALLOC_17(strlen(s)+1,char);
4106  strcpy(res,s);
4107  return res;
4108 }
4109 
static float swap_float(float source)
Definition: ioutils.cpp:206
MNE Surface or Volume (MneSurfaceOrVolume) class declaration.
#define FIFFV_MNE_COORD_MRI_VOXEL
MneProjData class declaration.
#define FIFF_MNE_SOURCE_SPACE_INTERPOLATOR
MneNearest class declaration.
IOUtils class declaration.
MneMghTagGroup class declaration.
FiffStream class declaration.
#define FIFF_BEM_SURF_TRIANGLES
Definition: fiff_file.h:739
#define FIFF_BEM_SURF_NNODE
Definition: fiff_file.h:736
Digitization points container and description.
One item in a derivation data set.
#define FIFF_MNE_SOURCE_SPACE_USE_TRIANGLES
MRI data volume geometry information like FreeSurfer keeps it.
Definition: mne_vol_geom.h:76
QSharedPointer< FiffDirNode > SPtr
Definition: fiff_dir_node.h:77
MNE Patch Information (MnePatchInfo) class declaration.
static bool fit_sphere_to_points(const Eigen::MatrixXf &rr, float simplex_size, Eigen::VectorXf &r0, float &R)
Digitization point description.
#define FIFF_MNE_SOURCE_SPACE_NUSE
#define FIFF_BEM_SURF_ID
Definition: fiff_file.h:734
Filter Thread Argument Description.
Sphere class declaration.
#define FIFF_MNE_SOURCE_SPACE_TYPE
This is used in the patch definitions.
Definition: mne_nearest.h:77
static qint32 swap_int(qint32 source)
Definition: ioutils.cpp:127
#define FIFF_MNE_SOURCE_SPACE_SELECTION
The MneMghTagGroup class.
#define FIFF_MNE_SOURCE_SPACE_TRIANGLES
static qint64 swap_long(qint64 source)
Definition: ioutils.cpp:161
#define FIFF_MNE_SOURCE_SPACE_NORMALS
#define FIFF_MNE_SOURCE_SPACE_NEAREST_DIST
#define FIFF_MNE_SOURCE_SPACE_NTRI
#define FIFF_MNE_SOURCE_SPACE_POINTS
QSharedPointer< FiffTag > SPtr
Definition: fiff_tag.h:152
Triangle data.
Definition: mne_triangle.h:75
MneMshDisplaySurface class declaration.
This defines a surface.
The MneMghTag class.
Definition: mne_mgh_tag.h:76
#define FIFF_MNE_SOURCE_SPACE_MRI_FILE
This defines a source space.
#define FIFFB_BEM
Definition: fiff_file.h:403
FIFF File I/O routines.
Definition: fiff_stream.h:104
Coordinate transformation descriptor.
#define FIFF_MNE_FILE_NAME
MneTriangle class declaration.
The MNE Msh Display Surface class holds information about a surface to be rendered.
#define FIFFB_BEM_SURF
Definition: fiff_file.h:404
#define FIFF_BEM_SURF_NTRI
Definition: fiff_file.h:737
#define FIFFV_MNE_COORD_RAS
#define FIFF_MNE_SOURCE_SPACE_VOXEL_DIMS
#define FIFF_MNE_SOURCE_SPACE_DIST
#define FIFF_BEM_SURF_NORMALS
Definition: fiff_file.h:740
FiffDigPoint class declaration.
QSharedPointer< FiffStream > SPtr
Definition: fiff_stream.h:107
#define FIFF_MNE_SOURCE_SPACE_DIST_LIMIT
fiff_float_t r[3]
FiffDigitizerData class declaration.
MneMghTag class declaration.
#define FIFF_MNE_SOURCE_SPACE_ID
#define FIFF_MNE_SOURCE_SPACE_NEAREST
Data associated with MNE computations for each mneMeasDataSet.
The MneProjData class.
Definition: mne_proj_data.h:74
static qint16 swap_short(qint16 source)
Definition: ioutils.cpp:114
MneSurfaceOld class declaration.
MneVolGeom class declaration.
MneSourceSpaceOld class declaration.
#define FIFF_SUBJ_HIS_ID
Definition: fiff_file.h:575
#define FIFF_MNE_SOURCE_SPACE_NPOINTS
Filter Thread Argument (FilterThreadArg) class declaration.
#define FIFF_BEM_COORD_FRAME
Definition: fiff_file.h:746
#define FIFF_BEM_SURF_NODES
Definition: fiff_file.h:738
#define FIFF_MNE_SOURCE_SPACE_NUSE_TRI
#define FIFF_BEM_SIGMA
Definition: fiff_file.h:747
#define FIFF_MNE_COORD_FRAME