[FV] = mesh_vertex_smooth(FV, index)
FV.vertices - Nx3
FV.faces - Nx3
FV.edges - NxN, see mesh_edges
index is the indices of FV.vertices to smooth (all by default)
This function shifts vertex location toward the center of mass of its
neighbouring vertices. It is useful to maintain surface smoothness.
This code is developed on the basis of Smith (2002, fig 5).
Smith, S. (2002). Fast robust automated brain extraction.
Human Brain Mapping, 17(3): 143-155.0001 function [FV] = mesh_vertex_smooth(FV, index), 0002 0003 % [FV] = mesh_vertex_smooth(FV, index) 0004 % 0005 % FV.vertices - Nx3 0006 % FV.faces - Nx3 0007 % FV.edges - NxN, see mesh_edges 0008 % 0009 % index is the indices of FV.vertices to smooth (all by default) 0010 % 0011 % This function shifts vertex location toward the center of mass of its 0012 % neighbouring vertices. It is useful to maintain surface smoothness. 0013 % 0014 % This code is developed on the basis of Smith (2002, fig 5). 0015 % 0016 % Smith, S. (2002). Fast robust automated brain extraction. 0017 % Human Brain Mapping, 17(3): 143-155. 0018 0019 0020 % This function adapts Smith, S. (2002), Fast robust automated brain 0021 % extraction. Human Brain Mapping, 17: 143-155. This function 0022 % corresponds to update component 2: surface smoothness control. 0023 0024 0025 if isfield(FV,'edges'), 0026 if isempty(FV.edges), 0027 FV.edges = mesh_edges(FV); 0028 end 0029 else 0030 FV.edges = mesh_edges(FV); 0031 end 0032 0033 % Define limits for radius of curvature, empirically optimized per 0034 % Smith (2002), see figure 6. 0035 radius_min = 3.33; % mm 0036 radius_max = 10.00; % mm 0037 % Sigmoid function parameters, 0038 % "where E and F control the scale and offset of the sigmoid" 0039 E = mean([(1 / radius_min), (1 / radius_max)]); 0040 F = 6 * ( (1 / radius_min) - (1 / radius_max) ); 0041 0042 Nvert = size(FV.vertices,1); 0043 0044 if ~exist('index','var'), index = 1:Nvert; end 0045 if isempty(index), index = 1:Nvert; end 0046 0047 fprintf('...smoothing %d vertices of surface triangulation\n',length(index)); 0048 0049 for i = 1:length(index), 0050 0051 vi = index(i); 0052 0053 v = FV.vertices(vi,:); 0054 0055 % Find neighbouring vertex coordinates 0056 neighbour_index = find(FV.edges(vi,:)); % the indices 0057 0058 % Find distances between vertex and neighbours, using edge lengths. 0059 % The mean value is lower(L) in Smith (2002, eq. 4) 0060 edge_distance = full(FV.edges(vi,neighbour_index)); 0061 L = mean(edge_distance); 0062 0063 % Find neighbour mean location; 0064 % this is 'mean position of A and B' 0065 % in figure 4 of Smith (2002) 0066 neighbour_vertices = FV.vertices(neighbour_index,:); 0067 Nmean = mean(neighbour_vertices); 0068 0069 % Find difference in distance between the central vertex and its 0070 % neighbours; this value is 's' in Smith (2002, fig 4); we take the 0071 % direction of this vector to be inward toward mean location of 0072 % neighbours 0073 s = Nmean - v; 0074 0075 % see function below 0076 vertNormal = surf_normal(FV, vi); 0077 0078 unit_normal = vertNormal / norm(vertNormal); 0079 0080 % Find the vector sn 0081 % the projection of s in the direction of the surface normal 0082 sn = dot( s, unit_normal ) * unit_normal; 0083 0084 % Find the vector st, the component of s orthogonal to the 0085 % surface normal vector. 0086 st = s - sn; 0087 0088 % Calculate radius of local curvature, solve Smith (2002, eq. 4) 0089 r = L^2 / (2 * norm(sn)); 0090 0091 % calculate the fraction of Sn to move the vertex 0092 f2 = (1 + tanh( F * (1 / r - E))) / 2; 0093 0094 % If we move the vertices at this point, all subsequent calculations 0095 % for the neighbours will account for the current move. So, we 0096 % calculate the surface normals during this loop, and we apply the 0097 % movements incrementally. If we accumulate the movement vectors and 0098 % apply them after the loop, the risk is that vertices will move too 0099 % far toward one another. 0100 0101 % sum([ v; st / 2; f2 * sn ]) 0102 0103 FV.vertices(vi,:) = v + (st / 2) + ( f2 * sn ); 0104 0105 plotface = 0; 0106 if plotface, 0107 figure; hold on 0108 %patch('faces',FV.faces,'vertices',FV.vertices,'facecolor',[.2 .2 .2],'facealpha',0.6) 0109 [faceIndexI,faceIndexJ] = find(FV.faces == vi); 0110 patch('faces',FV.faces(faceIndexI,:),'vertices',FV.vertices,'facecolor',[.6 .6 .6]) 0111 scatter3(v(1), v(2), v(3), 40, 'r', 'filled') 0112 tmp = neighbour_vertices; 0113 scatter3(tmp(:,1), tmp(:,2), tmp(:,3), 40, 'b', 'filled') 0114 vn = vertNormal; 0115 quiver3(v(1), v(2), v(3), vn(1), vn(2), vn(3), 0) 0116 tmpS = 10 * s; 0117 quiver3(v(1), v(2), v(3), tmpS(1), tmpS(2), tmpS(3), 0) 0118 tmpSn = 10 * sn; 0119 quiver3(v(1), v(2), v(3), tmpSn(1), tmpSn(2), tmpSn(3), 0) 0120 tmpSt = 10 * st; 0121 quiver3(v(1)+tmpSn(1), v(2)+tmpSn(2), v(3)+tmpSn(3), tmpSt(1), tmpSt(2), tmpSt(3), 0) 0122 0123 0124 motion = -10; 0125 tmpVn = v + (st / 2) + ( f2 * motion * sn); 0126 0127 0128 scatter3(tmpVn(1), tmpVn(2), tmpVn(3), 40, 'g', 'filled') 0129 rotate3d on 0130 end 0131 0132 end 0133 0134 % We could use St and Sn to move the vertex toward the center of mass 0135 % defined by the neighbour vertices, but this is best done inside the loop 0136 % above 0137 %F2 = repmat(f2,1,3); 0138 %FV.vertices = FV.vertices + (St / 2) + ( F2 .* Sn ); 0139 0140 return 0141 0142 0143 0144 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0145 function vertNormal = surf_normal(FV, vi), 0146 0147 % To calculate a vertex normal you need to perform an iteration 0148 % beginning with a triangle that holds the vertex and perform a cross 0149 % product on the two vectors of that triangle that extend from that 0150 % vertex. This will result in a normal for that triangle. Then go 0151 % around to each triangle containing that vertex and perform a cross 0152 % product for that triangle. You will end up with a set of normals for 0153 % all the triangles. To compute the vertex normal, sum all the face 0154 % normal vectors. 0155 0156 % get all the faces that contain this vertex 0157 [faceIndexI,faceIndexJ] = find(FV.faces == vi); 0158 0159 nface = length(faceIndexI); 0160 0161 faceNormals = zeros(nface,3); 0162 0163 for i = 1:nface, 0164 0165 f = faceIndexI(i); 0166 0167 i1 = FV.faces(f,1); 0168 i2 = FV.faces(f,2); 0169 i3 = FV.faces(f,3); 0170 0171 v1 = FV.vertices(i1,:); 0172 v2 = FV.vertices(i2,:); 0173 v3 = FV.vertices(i3,:); 0174 0175 edge1 = v2 - v1; 0176 edge2 = v3 - v1; 0177 0178 % If the vertices are given in clockwise order, when viewed from the 0179 % outside, the following calculates the "outward" surface normal. This 0180 % should be consistent with the matlab patch command. 0181 faceNormals(i,:) = cross( edge2, edge1 ); 0182 0183 % If the vertices are given in counterclockwise order, when viewed from 0184 % the outside, the following calculates the "outward" surface normal. 0185 % This should be consistent with the right hand rule. 0186 %faceNormals(i,:) = cross( edge1, edge2 ); 0187 0188 end 0189 0190 %faceNormalsMagnitude = vector_magnitude(faceNormals); 0191 %faceNormalsUnit = faceNormals ./ repmat(faceNormalsMagnitude,1,3); 0192 0193 % Area of Triangle = || AB x AC|| / 2 ; ie, the absolute value of 0194 % the length of the cross product of AB and AC divided by 2 0195 %faceArea = abs(faceNormalsMagnitude) / 2; 0196 0197 % Weight the faceNormals by the faceArea 0198 %vertNormals(v,:) = sum( faceNormals .* repmat(faceArea,1,3) ) / sum( faceArea ); 0199 0200 vertNormal = sum(faceNormals) / nface; 0201 0202 return