freesurfer_brain_correct - Correct intersections of brain/pial surfaces [FV_brain] = freesurfer_brain_correct(FV_brain,FV_pial) FV_brain.faces - Nx3 FV_brain.vertices - Nx3 FV_pial.faces - Nx3 FV_pial.vertices - Nx3 The FV_brain approximates an 'inner_skull', as it is the brain.tri file output by freesurfer during the skull strip procedure of the segmentation process. This surface can have intersections with the pial surface, which is disturbing, so we should take some trouble to shift it outward along the surface normal to avoid intersections. For practical purposes, we will have to do this regardless of the MRI intensity, so it is a detachment from the realistic geometry, unfortunatley. We might assume that a 2 mm separation is a reasonable minimum. The pial surface is assumed to be correct, so the 'brain' surface is moved outward along the surface normal, taking care not to introduce surface intersections in the brain surface itself (ie, constraining the inflation by neighbouring vertices). This is an iterative process for each vertex of the brain surface, each iteration moves a vertex and its neighbours by 0.5 mm.
0001 function [FV_brain] = freesurfer_brain_correct(FV_brain,FV_pial) 0002 0003 % freesurfer_brain_correct - Correct intersections of brain/pial surfaces 0004 % 0005 % [FV_brain] = freesurfer_brain_correct(FV_brain,FV_pial) 0006 % 0007 % FV_brain.faces - Nx3 0008 % FV_brain.vertices - Nx3 0009 % FV_pial.faces - Nx3 0010 % FV_pial.vertices - Nx3 0011 % 0012 % The FV_brain approximates an 'inner_skull', as it is the brain.tri 0013 % file output by freesurfer during the skull strip procedure of the 0014 % segmentation process. This surface can have intersections 0015 % with the pial surface, which is disturbing, so we should take 0016 % some trouble to shift it outward along the surface normal to avoid 0017 % intersections. For practical purposes, we will have to do this 0018 % regardless of the MRI intensity, so it is a detachment from the 0019 % realistic geometry, unfortunatley. We might assume that a 2 mm 0020 % separation is a reasonable minimum. The pial surface is assumed to 0021 % be correct, so the 'brain' surface is moved outward along the surface 0022 % normal, taking care not to introduce surface intersections in the brain 0023 % surface itself (ie, constraining the inflation by neighbouring vertices). 0024 % This is an iterative process for each vertex of the brain surface, each 0025 % iteration moves a vertex and its neighbours by 0.5 mm. 0026 % 0027 0028 % $Revision: 1.1 $ $Date: 2005/08/15 21:59:42 $ 0029 0030 % Copyright (C) 2005 Darren L. Weber 0031 % 0032 % This program is free software; you can redistribute it and/or 0033 % modify it under the terms of the GNU General Public License 0034 % as published by the Free Software Foundation; either version 2 0035 % of the License, or (at your option) any later version. 0036 % 0037 % This program is distributed in the hope that it will be useful, 0038 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0039 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0040 % GNU General Public License for more details. 0041 % 0042 % You should have received a copy of the GNU General Public License 0043 % along with this program; if not, write to the Free Software 0044 % Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, 0045 % USA. 0046 0047 % History: 09/2005, Darren.Weber_at_radiology.ucsf.edu 0048 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0049 0050 ver = '$Revision: 1.1 $ $Date: 2005/08/15 21:59:42 $'; 0051 fprintf('FREESURFER_BRAIN_CORRECT [v %s]\n',ver(11:15)); 0052 tic 0053 0054 %fprintf('...sorry, under development!\n\n'); return 0055 0056 DIST = 2; % inner skull to pial surface separation is 2 mm 0057 0058 DIST_MOVE = 0.25; % move 0.25 mm 0059 0060 FV = FV_brain; 0061 0062 if isfield(FV,'edges'), 0063 if isempty(FV.edges), 0064 FV.edges = mesh_edges(FV); 0065 end 0066 else 0067 FV.edges = mesh_edges(FV); 0068 end 0069 0070 0071 % get surface normals 0072 fprintf('...estimating brain vertex surface normals\n'); 0073 hf = figure('Visible','off'); 0074 hp = patch('faces',FV.faces,'vertices',FV.vertices); 0075 normals = get(hp,'VertexNormals'); 0076 close(hf); clear hf hp 0077 %Convert to unit normals 0078 [nrm,unit_normals] = colnorm(normals'); 0079 FV.normals = unit_normals'; 0080 clear nrm unit_normals 0081 0082 % % get surface normals 0083 % fprintf('...calculating vertex surface normals: pial\n'); 0084 % hf = figure('Visible','off'); 0085 % hp = patch('faces',FV_pial.faces,'vertices',FV_pial.vertices); 0086 % normals = get(hp,'VertexNormals'); 0087 % close(hf); clear hf hp 0088 % %Convert to unit normals 0089 % [nrm,normals] = colnorm(normals'); 0090 % FV_pial.normals = normals'; 0091 % clear nrm normals 0092 0093 indexSmallDist = ones(size(FV.vertices,1),1); 0094 0095 iteration = 0; 0096 while indexSmallDist, 0097 0098 iteration = iteration + 1; 0099 0100 % find vertices in the reduced tesselation that match those of the dense 0101 % tesselation; this can take several hours to run with dsearchn! 0102 fprintf('...find pial vertices nearest to each brain vertex\n'); 0103 indexPial4Brain = []; 0104 if exist('nearpoints','file'), 0105 [indexPial4Brain, bestDistSq] = nearpoints(FV.vertices',FV_pial.vertices'); 0106 indexPial4Brain = indexPial4Brain'; 0107 bestDist = sqrt(bestDistSq)'; 0108 else 0109 [indexPial4Brain, bestDist] = dsearchn(FV_pial.vertices,FV.vertices); 0110 bestDist = sqrt( bestDist .* bestDist ); 0111 end 0112 0113 0114 % find the brain vertices that are too close to the pial surface 0115 indexSmallDist = find(bestDist < DIST); 0116 0117 % find all the immediate neighbours of these vertices 0118 move_indices = []; 0119 for i = 1:length(indexSmallDist), 0120 index = indexSmallDist(i); 0121 nearest = find(FV.edges(:,index)); 0122 move_indices = [move_indices; index; nearest]; 0123 end 0124 move_indices = unique(move_indices); 0125 0126 % try to exclude the vertices of the cerebellum, 0127 % by excluding vertex distances > 6 mm 0128 tmpDist = bestDist(move_indices); 0129 tmpIndex = find(tmpDist > (3 * DIST) ); 0130 move_indices = setdiff(move_indices, tmpIndex); 0131 clear tmpDist tmpIndex 0132 0133 Nvert = length(move_indices); 0134 fprintf('...shifting %5d selected vertices outward %4.2f mm\n', Nvert, DIST_MOVE); 0135 0136 % We have an irregular surface (ie, the inner skull surface) and it 0137 % may not have an origin at (0,0,0). We have the surface vertex 0138 % normals, so now we move the vertex outward in the direction of 0139 % the vertex normal. We simply add a multiple of the unit vertex 0140 % normal to the current vertex xyz location. 0141 0142 move_vector = DIST_MOVE * FV.normals(move_indices,:); 0143 0144 move_vertices = FV.vertices(move_indices,:); 0145 0146 FV.vertices(move_indices,:) = move_vertices + move_vector; 0147 0148 if iteration > 250, 0149 warning('cannot correct brain surface'); 0150 return 0151 end 0152 0153 end 0154 0155 % apply one more move before smoothing 0156 move_vector = DIST_MOVE * FV.normals; 0157 move_vertices = FV.vertices; 0158 FV.vertices = move_vertices + move_vector; 0159 0160 % smooth the surface, which can shrink it slightly 0161 FV = mesh_vertex_smooth(FV); 0162 0163 FV_brain.vertices = FV.vertices; 0164 0165 t=toc; fprintf('...done (%5.2f sec).\n\n',t); 0166 0167 return 0168 0169 0170 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0171 % This code can be used to plot the vertices that are to be 0172 % corrected 0173 % 0174 % v = correct; 0175 % index = [scalpN oskullN]; 0176 % 0177 % patch('vertices',p.mesh.data.vertices{index(1)},'faces',p.mesh.data.faces{index(1)},... 0178 % 'FaceColor',[.9 .9 .9],'Edgecolor','none','FaceAlpha',.6); 0179 % daspect([1 1 1]); axis tight; hold on 0180 % 0181 % vert = p.mesh.data.vertices{index(1)}; 0182 % x = vert(v,1); 0183 % y = vert(v,2); 0184 % z = vert(v,3); 0185 % plot3(x,y,z,'ro') 0186 % 0187 % patch('vertices',p.mesh.data.vertices{index(2)},'faces',p.mesh.data.faces{index(2)},... 0188 % 'FaceColor',[.0 .6 .0],'Edgecolor',[.2 .2 .2],'FaceAlpha',.8); 0189 % 0190 % vert = p.mesh.data.vertices{index(2)}; 0191 % x = vert(v,1); 0192 % y = vert(v,2); 0193 % z = vert(v,3); 0194 % plot3(x,y,z,'bo')