mesh_vertex_smooth_check

 mesh_vertex_smooth_check(FV, index)

 FV.vertices - Nx3
 FV.faces - Nx3
 FV.edges - NxN, see mesh_edges

 index is any 1 vertex of FV.vertices (default is the first vertex)

 This function creates a graphic to check the vertex normal used in
 mesh_vertex_smooth. 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 mesh_vertex_smooth_check(FV, index),
0002 
0003 % mesh_vertex_smooth_check(FV, index)
0004 %
0005 % FV.vertices - Nx3
0006 % FV.faces - Nx3
0007 % FV.edges - NxN, see mesh_edges
0008 %
0009 % index is any 1 vertex of FV.vertices (default is the first vertex)
0010 %
0011 % This function creates a graphic to check the vertex normal used in
0012 % mesh_vertex_smooth. 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'),
0045     p = randperm(Nvert);
0046     index = p(1)
0047 end
0048 if isempty(index), 
0049     p = randperm(Nvert);
0050     index = p(1)
0051 end
0052 
0053 vi = index(1);
0054 
0055 v = FV.vertices(vi,:);
0056 
0057 % Find neighbouring vertex coordinates
0058 neighbour_index = find(FV.edges(vi,:));  % the indices
0059 
0060 % Find distances between vertex and neighbours, using edge lengths.
0061 % The mean value is lower(L) in Smith (2002, eq. 4)
0062 edge_distance = full(FV.edges(vi,neighbour_index));
0063 L = mean(edge_distance);
0064 
0065 % Find neighbour mean location;
0066 % this is 'mean position of A and B'
0067 % in figure 4 of Smith (2002)
0068 neighbour_vertices = FV.vertices(neighbour_index,:);
0069 Nmean = mean(neighbour_vertices);
0070 
0071 % Find difference in distance between the central vertex and its
0072 % neighbours; this value is 's' in Smith (2002, fig 4); we take the
0073 % direction of this vector to be inward toward mean location of
0074 % neighbours
0075 s = Nmean - v;
0076 
0077 % see function below
0078 vertNormal = surf_normal(FV, vi)
0079 
0080 unit_normal = vertNormal / norm(vertNormal);
0081 
0082 % Find the vector sn
0083 % the projection of s in the direction of the surface normal
0084 sn = dot( s, unit_normal ) * unit_normal;
0085 
0086 % Find the vector st, the component of s orthogonal to the
0087 % surface normal vector.
0088 st = s - sn;
0089 
0090 % Calculate radius of local curvature, solve Smith (2002, eq. 4)
0091 r = L^2 / (2 * norm(sn));
0092 
0093 % calculate the fraction of Sn to move the vertex
0094 f2 = (1 + tanh( F * (1 / r - E))) / 2;
0095 
0096 % If we move the vertices at this point, all subsequent calculations
0097 % for the neighbours will account for the current move.  So, we
0098 % calculate the surface normals during this loop, and we apply the
0099 % movements incrementally.  If we accumulate the movement vectors and
0100 % apply them after the loop, the risk is that vertices will move too
0101 % far toward one another.
0102 
0103 vn = sum([ v; st / 2; f2 * sn ])
0104 
0105 
0106 [faceIndexI,faceIndexJ] = find(FV.faces == vi);
0107 f = faceIndexI(1);
0108 
0109 vi1 = FV.faces(f,1);
0110 vi2 = FV.faces(f,2);
0111 vi3 = FV.faces(f,3);
0112 
0113 v1 = FV.vertices(vi1,:);
0114 v2 = FV.vertices(vi2,:);
0115 v3 = FV.vertices(vi3,:);
0116 
0117 vNorm = vertNormal;
0118 
0119 figure;
0120 subplot(1,2,1); hold on
0121 patch('faces',FV.faces,'vertices',FV.vertices,...
0122     'facecolor',[.8 .8 .8],'edgecolor',[.4 .4 .4],'facealpha',0.8)
0123 [faceIndexI,faceIndexJ] = find(FV.faces == vi);
0124 patch('faces',FV.faces(faceIndexI,:),'vertices',FV.vertices,'facecolor',[.6 .6 .6])
0125 scatter3(v(1), v(2), v(3), 60, 'r', 'filled')
0126 tmp = neighbour_vertices;
0127 scatter3(tmp(:,1), tmp(:,2), tmp(:,3), 40, 'b', 'filled')
0128 scatter3(v1(1), v1(2), v1(3), 40, 'r', 'filled')
0129 scatter3(v2(1), v2(2), v2(3), 40, 'g', 'filled')
0130 scatter3(v3(1), v3(2), v3(3), 40, 'y', 'filled')
0131 quiver3(v(1), v(2), v(3), vNorm(1), vNorm(2), vNorm(3),'r')
0132 rotate3d on
0133 
0134 
0135 % zoom in plot
0136 
0137 magnify = 10;
0138 vNorm = magnify * 10 * vertNormal;
0139 s = magnify * s;
0140 sn = magnify * sn;
0141 st = magnify * st;
0142 
0143 subplot(1,2,2); hold on
0144 scatter3(v(1), v(2), v(3), 40, 'r', 'filled')
0145 scatter3(vn(1), vn(2), vn(3), 40, 'g', 'filled')
0146 
0147 quiver3(v(1), v(2), v(3), vNorm(1), vNorm(2), vNorm(3),0,'r')
0148 quiver3(v(1), v(2), v(3), s(1), s(2), s(3),0,'g')
0149 quiver3(v(1), v(2), v(3), sn(1), sn(2), sn(3),0,'b')
0150 quiver3(v(1)+sn(1), v(2)+sn(2), v(3)+sn(3), st(1), st(2), st(3),0,'y')
0151 
0152 legend({'current vertex','new vertex','vertNorm','s','sn','st'},-1)
0153 
0154 [faceIndexI,faceIndexJ] = find(FV.faces == vi);
0155 patch('faces',FV.faces(faceIndexI,:),'vertices',FV.vertices,'facecolor',[.6 .6 .6])
0156 tmp = neighbour_vertices;
0157 scatter3(tmp(:,1), tmp(:,2), tmp(:,3), 40, 'b', 'filled')
0158 scatter3(v1(1), v1(2), v1(3), 40, 'r', 'filled')
0159 scatter3(v2(1), v2(2), v2(3), 40, 'c', 'filled')
0160 scatter3(v3(1), v3(2), v3(3), 40, 'y', 'filled')
0161 
0162 rotate3d on
0163 
0164 
0165 
0166 % We could use St and Sn to move the vertex toward the center of mass
0167 % defined by the neighbour vertices, but this is best done inside the loop
0168 % above
0169 %F2 = repmat(f2,1,3);
0170 %FV.vertices = FV.vertices + (St / 2) + ( F2 .* Sn );
0171 
0172 return
0173 
0174 
0175 
0176 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0177 function vertNormal = surf_normal(FV, vi),
0178 
0179 % To calculate a vertex normal you need to perform an iteration
0180 % beginning with a triangle that holds the vertex and perform a cross
0181 % product on the two vectors of that triangle that extend from that
0182 % vertex. This will result in a normal for that triangle. Then go
0183 % around to each triangle containing that vertex and perform a cross
0184 % product for that triangle. You will end up with a set of normals for
0185 % all the triangles. To compute the vertex normal, sum all the face
0186 % normal vectors.
0187 
0188 % get all the faces that contain this vertex
0189 [faceIndexI,faceIndexJ] = find(FV.faces == vi);
0190 
0191 nface = length(faceIndexI);
0192 
0193 faceNormals = zeros(nface,3);
0194 
0195 for i = 1:nface,
0196 
0197     f = faceIndexI(i);
0198 
0199     i1 = FV.faces(f,1);
0200     i2 = FV.faces(f,2);
0201     i3 = FV.faces(f,3);
0202 
0203     v1 = FV.vertices(i1,:);
0204     v2 = FV.vertices(i2,:);
0205     v3 = FV.vertices(i3,:);
0206 
0207     edge1 = v2 - v1;
0208     edge2 = v3 - v1;
0209     
0210     % If the vertices are given in clockwise order, when viewed from the
0211     % outside, the following calculates the "outward" surface normal. This
0212     % should be consistent with the matlab patch command.
0213     faceNormals(i,:) = cross( edge2, edge1 );
0214 
0215     % If the vertices are given in counterclockwise order, when viewed from
0216     % the outside, the following calculates the "outward" surface normal.
0217     % This should be consistent with the right hand rule.
0218     %faceNormals(i,:) = cross( edge1, edge2 );
0219 
0220 end
0221 
0222 %faceNormalsMagnitude = vector_magnitude(faceNormals);
0223 %faceNormalsUnit = faceNormals ./ repmat(faceNormalsMagnitude,1,3);
0224 
0225 % Area of Triangle = || AB x AC|| / 2 ; ie, the absolute value of
0226 % the length of the cross product of AB and AC divided by 2
0227 %faceArea = abs(faceNormalsMagnitude) / 2;
0228 
0229 % Weight the faceNormals by the faceArea
0230 %vertNormals(v,:) = sum( faceNormals .* repmat(faceArea,1,3) ) / sum( faceArea );
0231 
0232 vertNormal = sum(faceNormals) / nface;
0233 
0234 return