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_wfile0001 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