mesh_vertex_curvature

 This function adapts Smith, S. (2002), Fast robust automated brain
 extraction.  Human Brain Mapping, 17: 143-155.  This function
 corresponds to update component 2: surface smoothness control.

0001 function [FV] = mesh_vertex_smooth(FV,index,origin),
0002 
0003 % This function adapts Smith, S. (2002), Fast robust automated brain
0004 % extraction.  Human Brain Mapping, 17: 143-155.  This function
0005 % corresponds to update component 2: surface smoothness control.
0006 
0007 
0008 fprintf('this is in development\n');
0009 return
0010 
0011 
0012 xo = origin(1); yo = origin(2); zo = origin(3);
0013 
0014 v = FV.vertices(index,:);
0015 x = FV.vertices(index,1);
0016 y = FV.vertices(index,2);
0017 z = FV.vertices(index,3);
0018 
0019 % Find radial distance of vertex from origin
0020 r = sqrt( (x-xo)^2 + (y-yo)^2 + (z-zo)^2 );
0021 
0022 % Calculate unit vector
0023 v_unit_vector = ( v - origin ) / r;
0024 
0025 % Find direction cosines for line from center to vertex
0026 l = (x-xo)/r; % cos alpha
0027 m = (y-yo)/r; % cos beta
0028 n = (z-zo)/r; % cos gamma
0029 
0030 % Find neighbouring vertex coordinates
0031 vi = find(FV.edge(index,:));  % the indices
0032 neighbour_vertices = FV.vertices(vi,:);
0033 X = neighbour_vertices(:,1);
0034 Y = neighbour_vertices(:,2);
0035 Z = neighbour_vertices(:,3);
0036 
0037 % Find neighbour radial distances
0038 r_neighbours = sqrt( (X-xo).^2 + (Y-yo).^2 + (Z-zo).^2 );
0039 r_neighbours_mean = mean(r_neighbours);
0040 
0041 % Find difference in radial distance between the vertex of interest and its
0042 % neighbours; this value approximates the magnitude of sn in
0043 % Smith (2002, eq. 1 to 4)
0044 r_diff = r - r_neighbours_mean;
0045 
0046 % Find the vector sn, in the direction of the vertex of interest, given the
0047 % difference in radial distance between vertex and mean of neighbours
0048 sn = r_diff * v_unit_vector;
0049 
0050 % Find distances between vertex and neighbours, using edge lengths.
0051 % The mean value is l in Smith (2002, eq. 4)
0052 edge_distance = FV.edge(index,vi);
0053 edge_distance_mean = mean(edge_distance);
0054 
0055 % Calculate radius of local curvature, solve Smith (2002, eq. 4)
0056 if r_diff,
0057   radius_of_curvature = (edge_distance_mean ^ 2) / (2 * r_diff);
0058 else
0059   radius_of_curvature = 10000;
0060 end
0061 
0062 % Define limits for radius of curvature
0063 radius_min =  3.33; % mm
0064 radius_max = 10.00; % mm
0065 
0066 % Sigmoid function parameters,
0067 % "where E and F control the scale and offset of the sigmoid"
0068 E = mean([(1 / radius_min),  (1 / radius_max)]);
0069 F = 6 * ( (1 / radius_min) - (1 / radius_max) );
0070 
0071 Fsigmoid = (1 + tanh( F * (1 / radius_of_curvature - E))) / 2;
0072 
0073 % multiply sigmoid function by sn
0074 move_vector = Fsigmoid * sn;
0075 
0076 FV.vertices(index,:) = v + move_vector;
0077 
0078 return