Home > bioelectromagnetism > freesurfer_read_surf.m

freesurfer_read_surf

PURPOSE ^

freesurfer_read_surf - FreeSurfer I/O function to read a surface file

SYNOPSIS ^

function [vertices, faces] = freesurfer_read_surf(fname)

DESCRIPTION ^

 freesurfer_read_surf - FreeSurfer I/O function to read a surface file
 
 [vertices, faces] = freesurfer_read_surf(fname)
 
 Reads the vertex coordinates (mm) and face lists from a surface file.
 
 Surface files are stored as either triangulations or quadrangulations.
 That is, for a triangulation, each face is defined by 3 vertices.  For a
 quadrangulation, each face is defined by 4 vertices.  The rows of 'faces'
 contain indices into the rows of 'vertices', the latter holds the XYZ
 coordinates of each vertex.

 The freesurfer faces index the vertices in counter-clockwise order (when
 viewed from the outside of the surface).  This is consistent with a
 right-hand rule.  If we have vertices

 C           B


       A

 Then we can calculate an edge vector from A to B (ie, AB = B - A) and
 another edge vector from A to C (ie, AC = C - A).  If you form a "gun"
 with your thumb and forefinger of the right hand, then align your thumb
 with the AB vector and your forefinger with the AC vector, your palm is
 facing out of the screen and extending your middle finger in the
 orthogonal direction to the plane of the screen will give the outward
 surface normal of the triangle ABC.  (If you lookup "triangle" on
 Wolfram's mathworld, you can see that AB is referred to as c and AC is
 referred to as b.)

 However, if this surface is read into matlab, it will give INWARD surface
 normals in the matlab patch command.  For some reason, matlab is not
 following the right hand rule.  To get OUTWARD normals with the matlab
 patch command, use faces(:,[1 3 2]) (see below).

 The vertex coordinates are in mm.  The FreeSurfer coordinate
 system for surfaces is quite simple, but relating to their MRI
 cor-??? files is too confusing to explain here; see the FreeSurfer
 homepage or google the documentation by Graham Wideman.  For the
 surfaces, at least, the origin is somewhere in the center of the
 head, and the vertex XYZ coordinates are oriented such that +X is
 right, +Y is anterior and +Z is superior (this is the
 FreeSurfer RAS coordinate system).

 Note that reading the faces of a quad file can take a long
 time due to their compact storage format.  In this case, the return of
 vertices can be faster if the face output variable is not specified; in
 this case, the faces are not read.
 
 Try this to visualize the surface:
 Hp = patch('vertices',vertices,'faces',faces(:,[1 3 2]),...
       'facecolor',[.5 .5 .5],'edgecolor','none')
 camlight('headlight','infinite')
 vertnormals = get(Hp,'vertexnormals');

 See also freesurfer_write_surf, freesurfer_read_curv,
          freesurfer_read_wfile

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [vertices, faces] = freesurfer_read_surf(fname)
0002 
0003 % freesurfer_read_surf - FreeSurfer I/O function to read a surface file
0004 %
0005 % [vertices, faces] = freesurfer_read_surf(fname)
0006 %
0007 % Reads the vertex coordinates (mm) and face lists from a surface file.
0008 %
0009 % Surface files are stored as either triangulations or quadrangulations.
0010 % That is, for a triangulation, each face is defined by 3 vertices.  For a
0011 % quadrangulation, each face is defined by 4 vertices.  The rows of 'faces'
0012 % contain indices into the rows of 'vertices', the latter holds the XYZ
0013 % coordinates of each vertex.
0014 %
0015 % The freesurfer faces index the vertices in counter-clockwise order (when
0016 % viewed from the outside of the surface).  This is consistent with a
0017 % right-hand rule.  If we have vertices
0018 %
0019 % C           B
0020 %
0021 %
0022 %       A
0023 %
0024 % Then we can calculate an edge vector from A to B (ie, AB = B - A) and
0025 % another edge vector from A to C (ie, AC = C - A).  If you form a "gun"
0026 % with your thumb and forefinger of the right hand, then align your thumb
0027 % with the AB vector and your forefinger with the AC vector, your palm is
0028 % facing out of the screen and extending your middle finger in the
0029 % orthogonal direction to the plane of the screen will give the outward
0030 % surface normal of the triangle ABC.  (If you lookup "triangle" on
0031 % Wolfram's mathworld, you can see that AB is referred to as c and AC is
0032 % referred to as b.)
0033 %
0034 % However, if this surface is read into matlab, it will give INWARD surface
0035 % normals in the matlab patch command.  For some reason, matlab is not
0036 % following the right hand rule.  To get OUTWARD normals with the matlab
0037 % patch command, use faces(:,[1 3 2]) (see below).
0038 %
0039 % The vertex coordinates are in mm.  The FreeSurfer coordinate
0040 % system for surfaces is quite simple, but relating to their MRI
0041 % cor-??? files is too confusing to explain here; see the FreeSurfer
0042 % homepage or google the documentation by Graham Wideman.  For the
0043 % surfaces, at least, the origin is somewhere in the center of the
0044 % head, and the vertex XYZ coordinates are oriented such that +X is
0045 % right, +Y is anterior and +Z is superior (this is the
0046 % FreeSurfer RAS coordinate system).
0047 %
0048 % Note that reading the faces of a quad file can take a long
0049 % time due to their compact storage format.  In this case, the return of
0050 % vertices can be faster if the face output variable is not specified; in
0051 % this case, the faces are not read.
0052 %
0053 % Try this to visualize the surface:
0054 % Hp = patch('vertices',vertices,'faces',faces(:,[1 3 2]),...
0055 %       'facecolor',[.5 .5 .5],'edgecolor','none')
0056 % camlight('headlight','infinite')
0057 % vertnormals = get(Hp,'vertexnormals');
0058 %
0059 % See also freesurfer_write_surf, freesurfer_read_curv,
0060 %          freesurfer_read_wfile
0061 %
0062 
0063 % $Revision: 1.5 $ $Date: 2005/07/12 21:52:47 $
0064 
0065 % Copyright (C) 2000  Darren L. Weber
0066 %
0067 % This program is free software; you can redistribute it and/or
0068 % modify it under the terms of the GNU General Public License
0069 % as published by the Free Software Foundation; either version 2
0070 % of the License, or (at your option) any later version.
0071 %
0072 % This program is distributed in the hope that it will be useful,
0073 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0074 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0075 % GNU General Public License for more details.
0076 %
0077 % You should have received a copy of the GNU General Public License
0078 % along with this program; if not, write to the Free Software
0079 % Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,
0080 % USA.
0081 
0082 % History:  08/2000, Darren.Weber_at_radiology.ucsf.edu
0083 %
0084 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0085 
0086 ver = '$Revision: 1.5 $ $Date: 2005/07/12 21:52:47 $';
0087 fprintf('FREESURFER_READ_SURF [v %s]\n',ver(11:15));
0088 
0089 
0090 %QUAD_FILE_MAGIC_NUMBER =  (-1 & 0x00ffffff) ;
0091 %NEW_QUAD_FILE_MAGIC_NUMBER =  (-3 & 0x00ffffff) ;
0092 
0093 TRIANGLE_FILE_MAGIC_NUMBER  =  16777214;
0094 QUAD_FILE_MAGIC_NUMBER      =  16777215;
0095 
0096 
0097 % open it as a big-endian file
0098 fid = fopen(fname, 'rb', 'b');
0099 if (fid < 0),
0100     str = sprintf('could not open surface file %s.', fname);
0101     error(str);
0102 end
0103 
0104 fprintf('...reading surface file: %s\n', fname);
0105 tic;
0106 
0107 magic = freesurfer_fread3(fid);
0108 
0109 if (magic == QUAD_FILE_MAGIC_NUMBER),
0110     Nvertices = freesurfer_fread3(fid);
0111     Nfaces = freesurfer_fread3(fid);
0112     fprintf('...reading %d quad file vertices\n',Nvertices);
0113     vertices = fread(fid, Nvertices*3, 'int16') ./ 100 ; 
0114     if (nargout > 1),
0115         fprintf('...reading %d quad file faces (please wait)\n',Nfaces);
0116         faces = zeros(Nfaces,4);
0117         for iface = 1:Nfaces,
0118             for n=1:4,
0119                 faces(iface,n) = freesurfer_fread3(fid) ;
0120             end
0121             if(~rem(iface, 10000)), fprintf(' %7.0f',iface); end
0122             if(~rem(iface,100000)), fprintf('\n'); end
0123         end
0124     end
0125 elseif (magic == TRIANGLE_FILE_MAGIC_NUMBER),
0126     fprintf('...reading triangle file\n');
0127     tline = fgets(fid); % read creation date text line
0128     tline = fgets(fid); % read info text line
0129     
0130     Nvertices = fread(fid, 1, 'int32'); % number of vertices
0131     Nfaces = fread(fid, 1, 'int32'); % number of faces
0132     
0133     % vertices are read in column format and reshaped below
0134     vertices = fread(fid, Nvertices*3, 'float32');
0135     
0136     % faces are read in column format and reshaped
0137     faces = fread(fid, Nfaces*3, 'int32');
0138     faces = reshape(faces, 3, Nfaces)';
0139 else
0140     str = sprintf('unknown magic number in surface file %s.', fname);
0141     error(str);
0142 end
0143 
0144 vertices = reshape(vertices, 3, Nvertices)';
0145 fclose(fid);
0146 
0147 fprintf('...adding 1 to face indices for matlab compatibility.\n');
0148 faces = faces + 1;
0149 
0150 t=toc; fprintf('...done (%6.2f sec)\n\n',t);
0151 
0152 return

Generated on Mon 15-Aug-2005 15:36:19 by m2html © 2003