mne_read_bem_surfaces

 [surf] = mne_read_bem_surfaces(source,add_geom,tree)

 Reads source spaces from a fif file

 source      - The name of the file or an open file id
 add_geom    - Add geometry information to the surfaces
 tree        - Required if source is an open file id, search for the
               BEM surfaces here

0001 function [surf] = mne_read_bem_surfaces(source,add_geom,tree)
0002 %
0003 % [surf] = mne_read_bem_surfaces(source,add_geom,tree)
0004 %
0005 % Reads source spaces from a fif file
0006 %
0007 % source      - The name of the file or an open file id
0008 % add_geom    - Add geometry information to the surfaces
0009 % tree        - Required if source is an open file id, search for the
0010 %               BEM surfaces here
0011 %
0012 
0013 %
0014 %
0015 %   Author : Matti Hamalainen, MGH Martinos Center
0016 %   License : BSD 3-clause
0017 %
0018 %   Revision 1.3  2006/09/27 13:07:25  msh
0019 %   Added vertex normals to BEM surfaces.
0020 %   Fixed some data type casts in mne_read_surfaces
0021 %
0022 %   Revision 1.2  2006/05/03 18:53:06  msh
0023 %   Approaching Matlab 6.5 backward compatibility
0024 %
0025 %   Revision 1.1  2006/04/28 18:01:27  msh
0026 %   Added mne_read_bem_surfaces
0027 %   Fixed errors in mne_read_source_spaces and improved triangle data computation.
0028 %
0029 %
0030 
0031 me='MNE:mne_read_bem_surfaces';
0032 
0033 
0034 global FIFF;
0035 if isempty(FIFF)
0036    FIFF = fiff_define_constants();
0037 end
0038 
0039 if nargin == 1
0040     add_geom = false;
0041     fid      = -1;
0042     fname    = source;
0043 elseif nargin == 2
0044     fid      = -1;
0045     fname    = source;
0046 elseif nargin == 3
0047     a = whos('source');
0048     if strcmp(a.class,'char')
0049         fname = source;
0050         fid   = -1;
0051     else
0052         fid   = source;
0053     end
0054 else
0055     error(me,'Incorrect number of arguments');
0056 end
0057 %
0058 %   These fiff definitions are not needed elsewhere
0059 %
0060 FIFFB_BEM               = 310;    % BEM data
0061 FIFFB_BEM_SURF          = 311;    % One of the surfaces
0062 %
0063 FIFF_BEM_SURF_ID        = 3101;   % int    surface number
0064 FIFF_BEM_SURF_NAME      = 3102;   % string surface name
0065 FIFF_BEM_SURF_NNODE     = 3103;   % int    # of nodes on a surface
0066 FIFF_BEM_SURF_NTRI      = 3104;   % int    # number of triangles on a surface
0067 FIFF_BEM_SURF_NODES     = 3105;   % float  surface nodes (nnode,3)
0068 FIFF_BEM_SURF_TRIANGLES = 3106;   % int    surface triangles (ntri,3)
0069 FIFF_BEM_SURF_NORMALS   = 3107;   % float  surface node normal unit vectors (nnode,3)
0070 FIFF_BEM_COORD_FRAME    = 3112;   % The coordinate frame of the mode
0071 FIFF_BEM_SIGMA          = 3113;   % Conductivity of a compartment
0072 %
0073 %   Default coordinate frame
0074 %
0075 coord_frame = FIFF.FIFFV_COORD_MRI;
0076 %
0077 %   Open the file, create directory
0078 %
0079 if fid < 0
0080     [ fid, tree ] = fiff_open(fname);
0081     open_here = true;
0082 else
0083     open_here = false;
0084 end
0085 %
0086 %   Find BEM
0087 %
0088 bem = fiff_dir_tree_find(tree,FIFFB_BEM);
0089 if isempty(bem)
0090     close_file;
0091     error(me,'BEM data not found');
0092 end
0093 bem = bem(1);
0094 %
0095 %   Locate all surfaces
0096 %
0097 bemsurf = fiff_dir_tree_find(bem,FIFFB_BEM_SURF);
0098 if isempty(bemsurf)
0099     close_file;
0100     error(me,'BEM surface data not found');
0101 end
0102 fprintf(1,'\t%d BEM surfaces found\n',length(bemsurf));
0103 %
0104 %   Coordinate frame possibly at the top level
0105 %
0106 tag = find_tag(bem,FIFF_BEM_COORD_FRAME);
0107 if ~isempty(tag)
0108     coord_frame = tag.data;
0109 end
0110 %
0111 %   Read all surfaces
0112 %
0113 for k = 1:length(bemsurf)
0114     fprintf(1,'\tReading a surface...');
0115     this = read_source_space(bemsurf(k),coord_frame);
0116     fprintf(1,'[done]\n');
0117     if add_geom
0118         complete_surface_info();
0119     end
0120     surf(k) = this;
0121 end
0122 
0123 fprintf(1,'\t%d BEM surfaces read\n',length(surf));
0124 
0125 close_file();
0126 
0127 return;
0128 
0129     function [res] = read_source_space(this,def_coord_frame)
0130         %
0131         %   Read all the interesting stuff
0132         %
0133         tag = find_tag(this,FIFF_BEM_SURF_ID);
0134         if isempty(tag)
0135             res.id = FIFF.FIFFV_BEM_SURF_ID_UNKNOWN;
0136         else
0137             res.id = tag.data;
0138         end
0139 
0140         tag = find_tag(this,FIFF_BEM_SIGMA);
0141         if isempty(tag)
0142             res.sigma = 1.0;
0143         else
0144             res.sigma = tag.data;
0145         end
0146 
0147         tag = find_tag(this,FIFF_BEM_SURF_NNODE);
0148         if isempty(tag)
0149             close_file();
0150             error(me,'Number of vertices not found');
0151         end
0152         res.np = tag.data;
0153 
0154         tag = find_tag(this,FIFF_BEM_SURF_NTRI);
0155         if isempty(tag)
0156             close_file();
0157             error(me,'Number of triangles not found');
0158         else
0159             res.ntri = tag.data;
0160         end
0161 
0162         tag = find_tag(this,FIFF.FIFF_MNE_COORD_FRAME);
0163         if isempty(tag)
0164             tag = find_tag(this,FIFF_BEM_COORD_FRAME);
0165             if isempty(tag)
0166                 res.coord_frame = def_coord_frame;
0167             else
0168                 res.coord_frame = tag.data;
0169             end
0170         else
0171             res.coord_frame = tag.data;
0172         end
0173         %
0174         %   Vertices, normals, and triangles
0175         %
0176         tag = find_tag(this,FIFF_BEM_SURF_NODES);
0177         if isempty(tag)
0178             close_file();
0179             error(me,'Vertex data not found');
0180         end
0181         res.rr = tag.data;
0182         if size(res.rr,1) ~= res.np
0183             close_file();
0184             error(me,'Vertex information is incorrect');
0185         end
0186 
0187         tag = find_tag(this,FIFF.FIFF_MNE_SOURCE_SPACE_NORMALS);
0188         if isempty(tag)
0189             res.nn = [];
0190         else
0191             res.nn = tag.data;
0192             if size(res.nn,1) ~= res.np
0193                 close_file();
0194                 error(me,'Vertex normal information is incorrect');
0195             end
0196         end
0197         tag = find_tag(this,FIFF_BEM_SURF_TRIANGLES);
0198         if isempty(tag)
0199             close_file();
0200             error(me,'Triangulation not found');
0201         end
0202         res.tris = tag.data;
0203         if size(res.tris,1) ~= res.ntri
0204             close_file();
0205             error(me,'Triangulation information is incorrect');
0206         end
0207         return;
0208 
0209 
0210     end
0211 
0212     function [tag] = find_tag(node,findkind)
0213 
0214         for p = 1:node.nent
0215             if node.dir(p).kind == findkind
0216                 tag = fiff_read_tag(fid,node.dir(p).pos);
0217                 return;
0218             end
0219         end
0220         tag = [];
0221         return;
0222     end
0223 
0224     function complete_surface_info()
0225         %
0226         %   Main triangulation
0227         %
0228         fprintf(1,'\tCompleting triangulation info...');
0229         fprintf(1,'triangle normals...');
0230         this.tri_area = zeros(1,this.ntri);
0231         r1 = this.rr(this.tris(:,1),:);
0232         r2 = this.rr(this.tris(:,2),:);
0233         r3 = this.rr(this.tris(:,3),:);
0234         this.tri_cent = (r1+r2+r3)/3.0;
0235         this.tri_nn   = cross((r2-r1),(r3-r1));
0236         %
0237         %   Triangle normals and areas
0238         %
0239         for p = 1:this.ntri
0240             size = sqrt(this.tri_nn(p,:)*this.tri_nn(p,:)');
0241             this.tri_area(p) = size/2.0;
0242             if size > 0.0
0243                 this.tri_nn(p,:) = this.tri_nn(p,:)/size;
0244             end
0245         end
0246         %
0247         %   Accumulate the vertex normals
0248         %
0249         fprintf(1,'vertex normals...');
0250         this.nn = zeros(this.np,3);
0251         for p = 1:this.ntri
0252             this.nn(this.tris(p,:),:) = this.nn(this.tris(p,:),:) + kron(ones(3,1),this.tri_nn(p,:));
0253         end
0254         %
0255         %   Compute the lengths of the vertex normals and scale
0256         %
0257         fprintf(1,'normalize...');
0258         for p = 1:this.np
0259             size = sqrt(this.nn(p,:)*this.nn(p,:)');
0260             if size > 0
0261                 this.nn(p,:) = this.nn(p,:)/size;
0262             end
0263         end
0264         fprintf(1,'[done]\n');
0265     end
0266 
0267     function close_file()
0268 
0269         if open_here
0270             fclose(fid);
0271         end
0272 
0273     end
0274 
0275 
0276 end