emse_read_wfr

 emse_read_wfr - read EMSE wireframe (.wfr)

 [vertex,face,edge,meshtype] = emse_read_wfr(file,[options])

 All values returned are in meters.  With the returned structures, 
 create a vertex & face matrix:
 
   vertex_matrix = [vertex.x; vertex.y; vertex.z]';
   face_matrix = [face.vertex1; face.vertex2; face.vertex3]';

 These can be input to the patch command:

    Hpatch = patch('Vertices',vertex_matrix,'Faces',face_matrix,...
                   'EdgeColor',[.6 .6 .6],'FaceColor',[0.9 0.9 0.9]);

 See the patch and light commands to colour this object.
 
 'options' ... a cell array of strings.  By default it contains
 options = {'vertex','face','edge'}.  By default, this routine
 reads all available data from the emse file.  If 'options' is given
 and it doesn't contain one of these strings, that data will not be
 read or returned.

 meshtype is: 'unknown','scalp','outer skull','inner skull', or 'cortex'

0001 function [vertex,face,edge,meshtype] = emse_read_wfr(file,options)
0002 
0003 % emse_read_wfr - read EMSE wireframe (.wfr)
0004 %
0005 % [vertex,face,edge,meshtype] = emse_read_wfr(file,[options])
0006 %
0007 % All values returned are in meters.  With the returned structures,
0008 % create a vertex & face matrix:
0009 %
0010 %   vertex_matrix = [vertex.x; vertex.y; vertex.z]';
0011 %   face_matrix = [face.vertex1; face.vertex2; face.vertex3]';
0012 %
0013 % These can be input to the patch command:
0014 %
0015 %    Hpatch = patch('Vertices',vertex_matrix,'Faces',face_matrix,...
0016 %                   'EdgeColor',[.6 .6 .6],'FaceColor',[0.9 0.9 0.9]);
0017 %
0018 % See the patch and light commands to colour this object.
0019 %
0020 % 'options' ... a cell array of strings.  By default it contains
0021 % options = {'vertex','face','edge'}.  By default, this routine
0022 % reads all available data from the emse file.  If 'options' is given
0023 % and it doesn't contain one of these strings, that data will not be
0024 % read or returned.
0025 %
0026 % meshtype is: 'unknown','scalp','outer skull','inner skull', or 'cortex'
0027 %
0028 
0029 % $Revision: 1.1 $ $Date: 2005/01/21 04:33:18 $
0030 
0031 % Copyright (C) 2001  Darren L. Weber
0032 %
0033 % This program is free software; you can redistribute it and/or
0034 % modify it under the terms of the GNU General Public License
0035 % as published by the Free Software Foundation; either version 2
0036 % of the License, or (at your option) any later version.
0037 %
0038 % This program is distributed in the hope that it will be useful,
0039 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0040 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0041 % GNU General Public License for more details.
0042 %
0043 % You should have received a copy of the GNU General Public License
0044 % along with this program; if not, write to the Free Software
0045 % Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,
0046 % USA.
0047 
0048 % History:  10/98 Abbas Kouzani (egazk@flinders.edu.au)
0049 %           09/01 Darren.Weber_at_radiology.ucsf.edu
0050 %                 - created function, rather than script
0051 %                 - added functionality to handle different
0052 %                   minor revisions.
0053 %           03/02 Darren.Weber_at_radiology.ucsf.edu
0054 %                 - optimised fscanf processes for matlab, the
0055 %                   function now operates in less than 1/2 time,
0056 %                   especially for revision 3 data.
0057 %           12/04 Darren.Weber_at_radiology.ucsf.edu
0058 %                 - modified function name from mesh_emse2matlab
0059 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0060 
0061 ver = '$Revision: 1.1 $ $Date: 2005/01/21 04:33:18 $';
0062 fprintf('\nEMSE_READ_WFR [v%s]\n',ver(11:15));
0063 
0064 if ~exist('options','var'),
0065     options = {'vertex','face','edge'};
0066 end
0067 
0068 [path,name,ext] = fileparts(file);
0069 file = fullfile(path,[name ext]);
0070 
0071 [fid,msg] = fopen(file,'r');
0072 if ~isempty(msg), error(msg); end
0073 
0074 tic;
0075 
0076 % Read prolog
0077 version   =fscanf(fid,'%f',1);
0078 file_type =fscanf(fid,'%f',1);
0079 minor_rev =fscanf(fid,'%f',1);
0080 
0081 fprintf('...WFR version = %d, minor_revision = %d, file-type = %d\n',...
0082     version,minor_rev,file_type);
0083 
0084 if ~(file_type==8000 | file_type==4000)
0085     msg = sprintf('cannot read WFR file type: %d',file_type);
0086     error(msg);
0087 end
0088 
0089 % Read header (format depends on minor revision)
0090 if(minor_rev==3)
0091     type      =fscanf(fid,'%f',1);
0092 else
0093     radius    =fscanf(fid,'%f',1);
0094     vertex_no =fscanf(fid,'%d',1);
0095     face_no   =fscanf(fid,'%d',1);
0096     edge_no   =fscanf(fid,'%d',1);
0097     
0098     fprintf('...radius = %f meters\n',radius);
0099     fprintf('...reading %d vertices\n...%d faces\n...%d edges\n',...
0100         vertex_no,face_no,edge_no);
0101 
0102     if(minor_rev==1), type = 0;
0103     else,             type = fscanf(fid,'%f',1);
0104     end
0105 end
0106 
0107 if(minor_rev==1)
0108     meshtype = 'unknown';
0109 else
0110     switch type
0111         case 0,             meshtype = 'unknown';
0112         case { 64,  40},    meshtype = 'scalp';
0113         case {128,  80},    meshtype = 'outer skull';
0114         case {256, 100},    meshtype = 'inner skull';
0115         case {512, 200},    meshtype = 'cortex';
0116         otherwise,          meshtype = 'unknown';
0117     end
0118 end
0119 
0120 fprintf('...mesh type: %s\n',meshtype);
0121 
0122 % Read data (format depends on minor revision)
0123 vertex = struct([]);
0124 face = struct([]);
0125 edge = struct([]);
0126 
0127 if(minor_rev==3)
0128     
0129     % Read the whole file
0130     fprintf('...reading minor revision 3 data');
0131     Tmp = fscanf(fid,'%s%f%f%f',[4,inf]);
0132     fprintf('...done\n');
0133     
0134     % Vertices
0135     if(max(strcmp('vertex',options)) > 0),
0136         fprintf('...creating vertex struct');
0137         % first get the numeric code for 'v'
0138         vcode = Tmp(1,1);
0139         % now find all columns of Tmp with vcode in 1st row
0140         vindex = find(Tmp(1,:) == vcode);
0141         
0142         tmp = Tmp(2:4,vindex);
0143         tmp = num2cell(tmp);
0144         vertex = struct(...
0145             'index', [1:length(vindex)],...
0146             'x',     tmp(1,:),...
0147             'y',     tmp(2,:),...
0148             'z',     tmp(3,:));
0149         clear tmp;
0150         fprintf('...done\n');
0151     else
0152         fprintf('...skipping vertices\n');
0153     end
0154     
0155     % Faces
0156     if(max(strcmp('face',options)) > 0),
0157         fprintf('...creating face struct');
0158         % first get the numeric code for 't'
0159         tcode = Tmp(1,length(vindex)+1);
0160         % now find all columns of Tmp with tcode in 1st row
0161         tindex = find(Tmp(1,:) == tcode);
0162         
0163         % matlab vertex indices start at one,
0164         % not zero, so we add one to these emse values
0165         tmp = Tmp(2:4,tindex) + 1;
0166         
0167         tmp = num2cell(tmp);
0168         face = struct(...
0169             'index',   [1:length(tindex)],...
0170             'vertex1', tmp(1,:),...
0171             'vertex2', tmp(2,:),...
0172             'vertex3', tmp(3,:));
0173         clear tmp;
0174         fprintf('...done\n');
0175     else
0176         fprintf('...skipping faces\n');
0177     end
0178     
0179     % Edges
0180     fprintf('...no edges for minor revision 3\n');
0181     clear Tmp;
0182     
0183 elseif(minor_rev~=4)
0184     
0185     % minor revision 1 & 2 format
0186     disp('...reading minor revision 1 or 2 data');
0187     if(max(strcmp('vertex',options)) > 0),
0188         fprintf('...reading %d vertices',vertex_no);
0189         
0190         tmp = zeros(13,vertex_no);
0191         tmp = fscanf(fid,'%d%x%d%d%g%g%g%g%g%g%g%g%g',[13,vertex_no]);
0192         
0193         tmp(1,:) = tmp(1,:) + 1;
0194         tmp      = num2cell(tmp);
0195         vertex   = struct(...
0196             'index',        tmp( 1,:),...
0197             'address',      tmp( 2,:),...
0198             'channel_index',tmp( 3,:),...
0199             'x',            tmp( 5,:),...
0200             'y',            tmp( 6,:),...
0201             'z',            tmp( 7,:),...
0202             'xnormal',      tmp( 9,:),...
0203             'ynormal',      tmp(10,:),...
0204             'znormal',      tmp(11,:),...
0205             'potential',    tmp(12,:),...
0206             'curvature',    tmp(13,:));
0207         clear tmp;
0208         fprintf('...done\n');
0209     else
0210         fprintf('...skipping vertices\n');
0211         tmp = zeros(13,vertex_no);
0212         tmp = fscanf(fid,'%d%x%d%d%g%g%g%g%g%g%g%g%g',[13,vertex_no]);
0213         clear tmp;
0214     end
0215     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0216     if(max(strcmp('face',options)) > 0),
0217         fprintf('...reading %d faces',face_no);
0218         
0219         tmp = zeros(18,face_no);
0220         tmp = fscanf(fid,'%d%x%g%g%g%g%g%g%g%g%g%g%x%x%x%x%x%x',[18,face_no]);
0221         tmp(1,:) = tmp(1,:) + 1;
0222         tmp      = num2cell(tmp);
0223         face     = struct(...
0224             'index',        tmp( 1,:),...
0225             'address',      tmp( 2,:),...
0226             'solid_angle',  tmp( 3,:),...
0227             'magnitude',    tmp( 4,:),...
0228             'potential',    tmp( 5,:),...
0229             'area',         tmp( 6,:),...
0230             'center_x',     tmp( 7,:),...
0231             'center_y',     tmp( 8,:),...
0232             'center_z',     tmp( 9,:),...
0233             'normal_x',     tmp(10,:),...
0234             'normal_y',     tmp(11,:),...
0235             'normal_z',     tmp(12,:),...
0236             'vertex1',      tmp(13,:),...
0237             'vertex2',      tmp(14,:),...
0238             'vertex3',      tmp(15,:),...
0239             'edge1',        tmp(16,:),...
0240             'edge2',        tmp(17,:),...
0241             'edge3',        tmp(18,:));
0242         clear tmp;
0243         fprintf('...done\n');
0244         
0245         % In minor rev4, the face vertex and edges
0246         % refer to the vertex.address field, so this
0247         % is corrected here.  Not sure how to avoid 'for'
0248         fprintf('...converting face vertices from address to index (this takes a while)');
0249         for i=1:face_no,
0250             face(i).vertex1 = find([vertex.address] == face(i).vertex1);
0251             face(i).vertex2 = find([vertex.address] == face(i).vertex2);
0252             face(i).vertex3 = find([vertex.address] == face(i).vertex3);
0253         end
0254         fprintf('...done\n');
0255         
0256         if(max(strcmp('edge',options)) > 0),
0257             fprintf('...converting edges from address to index (this takes a while)');
0258             for i=1:face_no,
0259                 face(i).edge1 = find([vertex.address] == face(i).edge1);
0260                 face(i).edge2 = find([vertex.address] == face(i).edge2);
0261                 face(i).edge3 = find([vertex.address] == face(i).edge3);
0262             end
0263             fprintf('...done\n');
0264         end
0265     else
0266         disp('...skipping faces');
0267         tmp = zeros(18,face_no);
0268         tmp = fscanf(fid,'%d%x%g%g%g%g%g%g%g%g%g%g%x%x%x%x%x%x',[18,face_no]);
0269         clear tmp;
0270     end
0271     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0272     if(max(strcmp('edge',options)) > 0),
0273         fprintf('...reading %d edges',edge_no);
0274         
0275         tmp = zeros(4,edge_no);
0276         tmp = fscanf(fid,'%d%x%x%x',[4,edge_no]);
0277         
0278         tmp(1,:) = tmp(1,:) + 1;
0279         tmp      = num2cell(tmp);
0280         edge     = struct(...
0281             'index',        tmp(1,:),...
0282             'address',      tmp(2,:),...
0283             'vertex1',      tmp(3,:),...
0284             'vertex2',      tmp(4,:));
0285         clear tmp;
0286         fprintf('...done\n');
0287         fprintf('...converting edge vertices from address to index (this takes a while)');
0288         for i=1:edge_no,
0289             edge(i).vertex1 = find([vertex.address] == edge(i).vertex1);
0290             edge(i).vertex2 = find([vertex.address] == edge(i).vertex2);
0291         end
0292         fprintf('...done\n');
0293     else
0294         disp('...skipping edges');
0295         %for i=1:edge_no,
0296         %    tmp=fscanf(fid,'%*d',1);
0297         %    tmp=fscanf(fid,'%*x',3);
0298         %end
0299     end
0300 else
0301     % minor revision 4 format
0302     disp('...Reading Minor Revision 4 Data');
0303     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0304     if(max(strcmp('vertex',options)) > 0),
0305         fprintf('...reading %d vertices',vertex_no);
0306         
0307         index = meshgrid(1:1:vertex_no,1);
0308         index = num2cell(index);
0309         
0310         tmp = zeros(11,vertex_no);
0311         tmp = fscanf(fid,'%d%d%g%g%g%d%g%g%g%g%g',[11,vertex_no]);
0312         
0313         tmp      = num2cell(tmp);
0314         vertex   = struct(...
0315             'index',        index,...
0316             'channel_index',tmp( 1,:),...
0317             'x',            tmp( 3,:),...
0318             'y',            tmp( 4,:),...
0319             'z',            tmp( 5,:),...
0320             'xnormal',      tmp( 7,:),...
0321             'ynormal',      tmp( 8,:),...
0322             'znormal',      tmp( 9,:),...
0323             'potential',    tmp(10,:),...
0324             'curvature',    tmp(11,:));
0325         clear tmp index;
0326         fprintf('...done\n');
0327     else
0328         disp('...skipping vertices');
0329         tmp = zeros(11,vertex_no);
0330         tmp = fscanf(fid,'%d%d%g%g%g%d%g%g%g%g%g',[11,vertex_no]);
0331     end
0332     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0333     if(max(strcmp('face',options)) > 0),
0334         fprintf('...reading %d faces',face_no);
0335         
0336         index = meshgrid(1:1:face_no,1);
0337         index = num2cell(index);
0338         
0339         tmp = zeros(16,face_no);
0340         tmp = fscanf(fid,'%g%g%g%g%g%g%g%g%g%g%d%d%d%d%d%d',[16,face_no]);
0341         
0342         tmp(11:16,:) = tmp(11:16,:) + 1;
0343         
0344         tmp      = num2cell(tmp);
0345         face     = struct(...
0346             'index',        index,...
0347             'solid_angle',  tmp( 1,:),...
0348             'magnitude',    tmp( 2,:),...
0349             'potential',    tmp( 3,:),...
0350             'area',         tmp( 4,:),...
0351             'center_x',     tmp( 5,:),...
0352             'center_y',     tmp( 6,:),...
0353             'center_z',     tmp( 7,:),...
0354             'normal_x',     tmp( 8,:),...
0355             'normal_y',     tmp( 9,:),...
0356             'normal_z',     tmp(10,:),...
0357             'vertex1',      tmp(11,:),...
0358             'vertex2',      tmp(12,:),...
0359             'vertex3',      tmp(13,:),...
0360             'edge1',        tmp(14,:),...
0361             'edge2',        tmp(15,:),...
0362             'edge3',        tmp(16,:));
0363         clear tmp index;
0364         fprintf('...done\n');
0365     else
0366         disp('...skipping faces');
0367         tmp = zeros(16,face_no);
0368         tmp = fscanf(fid,'%g%g%g%g%g%g%g%g%g%g%d%d%d%d%d%d',[16,face_no]);
0369     end
0370     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0371     if(max(strcmp('edge',options)) > 0),
0372         fprintf('...reading %d edges',edge_no);
0373         
0374         index = meshgrid(1:1:edge_no,1);
0375         address = num2cell(index * 0);
0376         index = num2cell(index);
0377         
0378         tmp = zeros(2,edge_no);
0379         tmp = fscanf(fid,'%d%d',[2,edge_no]);
0380         tmp = tmp + 1;
0381         tmp = num2cell(tmp);
0382         
0383         edge = struct(...
0384             'index',   index,...
0385             'address', address,...
0386             'vertex1', tmp(1,:),...
0387             'vertex2', tmp(2,:));
0388         clear tmp index address;
0389         fprintf('...done\n');
0390     else
0391         fprintf('...skipping edges\n');
0392         %tmp=fscanf(fid,'%*d',2*edge_no);
0393     end
0394 end
0395 fclose(fid);
0396 
0397 t=toc; fprintf('...done (%5.2f sec)\n\n',t);
0398 
0399 return