avw_bet

 avw_bet

 [ u ] = avw_bet(avw,bt)

 avw - a struct returned by avw_read
 bt - this modulates the intensity weighting of the fit, values vary from
 0-1, the default is 0.5 (larger values give ?larger? brain).

 This code is developed entirely on the basis of Smith, S. (2002). Fast
 robust automated brain extraction. Human Brain Mapping, 17(3): 143-155.

0001 function [ FV ] = avw_bet(avw,bt)
0002 
0003 % avw_bet
0004 %
0005 % [ u ] = avw_bet(avw,bt)
0006 %
0007 % avw - a struct returned by avw_read
0008 % bt - this modulates the intensity weighting of the fit, values vary from
0009 % 0-1, the default is 0.5 (larger values give ?larger? brain).
0010 %
0011 % This code is developed entirely on the basis of Smith, S. (2002). Fast
0012 % robust automated brain extraction. Human Brain Mapping, 17(3): 143-155.
0013 %
0014 
0015 
0016 % see Smith (2002, p. 150)
0017 % bt can vary between 0-1
0018 if ~exist('bt','var'), bt = 0.5; end
0019 if isempty(bt), bt = 0.5; end
0020 
0021 fprintf('...estimating brain surface...'); tic;
0022 
0023 [FV,t,t02,t98,tm,COG] = avw_bet_init(avw);
0024 
0025 Hf = figure;
0026 Hp = patch('vertices',FV.vertices,...
0027     'faces',FV.faces,...
0028     'facecolor',[.8 .75 .75],...
0029     'edgecolor',[.6 .6 .6]);
0030 light
0031 drawnow
0032 
0033 for i = 1:1000,
0034     
0035     updateStr = sprintf('update iteration: %04d\n',i);
0036     set(Hf,'Name',['avw_bet, ',updateStr]);
0037     drawnow
0038 
0039     % calculate the surface normals
0040     [FV.normals,FV.unit_normals] = mesh_vertex_normals(FV);
0041 
0042     % component 1
0043     [S,Sn,St,Sn_unit] = avw_bet_update1(FV);
0044 
0045     % component 2
0046     [f2,L] = avw_bet_update2(FV,COG.voxels,Sn);
0047 
0048     % component 3
0049     [ f3, t1, Imin, Imax] = avw_bet_update3(FV,avw,bt,Sn,t,t02,tm);
0050 
0051 
0052     % Nvert is size(FV.vertices,1)
0053     % st is a vector Nvert x 3 (tangential to local surface)
0054     % sn is a vector Nvert x 3 (normal to local surface)
0055     % sn_unit is the unit vector of sn, a vector Nvert x 3
0056     % f3 is a scalar Nvert x 1 (fractional update of intensity)
0057     % L  is a scalar Nvert x 1 (the mean distance of a vertex to its neighbours)
0058 
0059 
0060     % check on the units of sn and L - voxels or mm?
0061 
0062     f2 = repmat(f2,1,3);
0063     f3 = repmat(f3,1,3);
0064     L  = repmat(L ,1,3);
0065 
0066     u1 = 0.5 * St;
0067     u2 = f2 .* Sn;
0068     u3 = 0.05 .* f3 .* L .* Sn_unit;
0069 
0070     u = u1 + u2 + u3;
0071 
0072     % Now apply movements to FV.vertices - somehow?
0073     FV.vertices = FV.vertices + u;
0074 
0075     set(Hp,'vertices',FV.vertices);
0076     drawnow
0077     
0078 end
0079 
0080 
0081 % Now mask avw.img to remove all voxels outside of FV
0082 % This will be tricky.
0083 
0084 
0085 
0086 
0087 time = toc; fprintf('done (%5.2f sec)\n',time);
0088 return
0089 
0090 
0091 
0092 
0093 
0094 
0095 %------------------------------------------------------
0096 function [FV,t,t02,t98,tm,COG] = avw_bet_init(avw)
0097 
0098 % t02 and t98 are 2% and 98% of cumuative volume histogram
0099 
0100 [bins,freq,freq_nozero] = avw_histogram(avw,2,[],0);
0101 
0102 
0103 % this seems to work, but not sure if it is a correct definition of
0104 % 'cumulative volume histogram'
0105 binbyfreq = bins .* freq;
0106 
0107 t02_tmp = 0.02 * sum(binbyfreq);
0108 t98_tmp = 0.98 * sum(binbyfreq);
0109 
0110 i = 1;
0111 while sum(binbyfreq(1:i)) < t02_tmp, i = i + 1; end
0112 t02 = bins(i-1);
0113 
0114 i = 1; while sum(binbyfreq(1:i)) < t98_tmp, i = i + 1; end
0115 t98 = bins(i);
0116 
0117 % t is the 'brain/background' threshold, which is used to roughly estimate
0118 % the position of the center of gravity of the brain/head in the image
0119 % volume
0120 
0121 % t is simply set to lie 10% of the way between t02 and t98
0122 t = (0.1 * (t98 - t02)) + t02;
0123 
0124 %Create binarized image to determine Center of gravity (COG) and radius of
0125 %sphere volume to determine tm.
0126 bin = avw;
0127 % find all voxels with intensity less than t
0128 voxels_lt_t = find(avw.img < t);
0129 bin.img(voxels_lt_t) = 0;
0130 % find all voxels with intensity greater than t98
0131 voxels_gt_t98 = find(avw.img > t98);
0132 bin.img(voxels_gt_t98) = 0;
0133 % binary image
0134 nonzero = find(bin.img);
0135 bin.img(nonzero) = 1;
0136 
0137 % find all voxel indices for these thresholded voxels
0138 [i,j,k] = ind2sub(size(avw.img),nonzero);
0139 
0140 
0141 % find weighted sum of positions (CENTER OF GRAVITY (COG))
0142 COG = avw_center_mass(bin);
0143 xCOG = COG.voxels(1);
0144 yCOG = COG.voxels(2);
0145 zCOG = COG.voxels(3);
0146 
0147 radius_voxel = sqrt( (i-xCOG).^2 + (j-yCOG).^2 + (k-zCOG).^2 );
0148 
0149 % xCOG = COG.mm(1);
0150 % yCOG = COG.mm(2);
0151 % zCOG = COG.mm(3);
0152 %
0153 % i = i * double(avw.hdr.dime.pixdim(2));
0154 % j = j * double(avw.hdr.dime.pixdim(3));
0155 % k = k * double(avw.hdr.dime.pixdim(4));
0156 %
0157 % radius_mm = sqrt( (i-xCOG).^2 + (j-yCOG).^2 + (k-zCOG).^2 );
0158 
0159 radius_mean = mean(radius_voxel) / 2;
0160 
0161 % "Finally, the median intensity of all points within a sphere of the
0162 % estimated radius and centered on the estimated COG is found (tm)"
0163 % Smith(2002, p. 146).
0164 
0165 tm = median(avw.img(nonzero));
0166 
0167 % Initialise the spherical tesselation
0168 % with 4 recursions we get 2562 vertices,
0169 % with 5 recursions we get 10,242 vertices;
0170 FV = sphere_tri('ico',4,radius_mean);
0171 FV.edge = mesh_edges(FV);
0172 
0173 Nvert = size(FV.vertices,1);
0174 A = repmat([xCOG,yCOG,zCOG],Nvert,1);
0175 FV.vertices = FV.vertices + A;
0176 
0177 return
0178 
0179 
0180 
0181 
0182 
0183 
0184 %------------------------------------------------------
0185 function [S,Sn,St,Sn_unit] = avw_bet_update1(FV),
0186 
0187 % This function implements Smith, S. (2002), Fast robust automated brain
0188 % extraction.  Human Brain Mapping, 17: 143-155.  It corresponds to update
0189 % component 1: within surface vertex spacing.
0190 
0191 Nvert = size(FV.vertices,1);
0192 
0193 for index = 1:Nvert,
0194 
0195     v = FV.vertices(index,:);
0196     x = v(1);
0197     y = v(2);
0198     z = v(3);
0199 
0200     unit_normal = FV.unit_normals(index,:);
0201 
0202     % Find neighbouring vertex coordinates
0203     vi = find(FV.edge(index,:));  % the indices
0204     neighbour_vertices = FV.vertices(vi,:);
0205     X = neighbour_vertices(:,1);
0206     Y = neighbour_vertices(:,2);
0207     Z = neighbour_vertices(:,3);
0208 
0209     % Find neighbour mean location; this is 'mean position of A and B' in
0210     % figure 4 of Smith (2002)
0211     Xmean = mean(X);
0212     Ymean = mean(Y);
0213     Zmean = mean(Z);
0214 
0215     % Find difference in distance between the vertex of interest and its
0216     % neighbours; this value is 's' and 'sn' in figure 4 of
0217     % Smith (2002, eq. 1 to 4)
0218     s = [ Xmean - x, Ymean - y, Zmean - z]; % inward toward mean
0219 
0220     % Find the vector sn
0221     sn = dot( s, unit_normal ) * unit_normal;
0222 
0223     % Find the vector st
0224     st =  s - sn; % absolute value
0225 
0226     S(index,:) = s;
0227     Sn(index,:) = sn;
0228     St(index,:) = st;
0229 
0230     Sn_unit(index,:) = sn ./ sqrt( (sn(1) - v(1)).^2 + (sn(2) - v(2)).^2 + (sn(3) - v(3)).^2 );
0231     
0232 end
0233 
0234 return
0235 
0236 
0237 
0238 
0239 
0240 
0241 %------------------------------------------------------
0242 function [f2,L] = avw_bet_update2(FV,origin,Sn),
0243 
0244 % This function adapts Smith, S. (2002), Fast robust automated brain
0245 % extraction.  Human Brain Mapping, 17: 143-155.  This function
0246 % corresponds to update component 2: surface smoothness control.
0247 
0248 xo = origin(1); yo = origin(2); zo = origin(3);
0249 
0250 % Define limits for radius of curvature empirically optimized per Smith
0251 % (2002), see figure 6.
0252 radius_min =  3.33; % mm
0253 radius_max = 10.00; % mm
0254 % Sigmoid function parameters,
0255 % "where E and F control the scale and offset of the sigmoid"
0256 E = mean([(1 / radius_min),  (1 / radius_max)]);
0257 F = 6 * ( (1 / radius_min) - (1 / radius_max) );
0258 
0259 Nvert = size(FV.vertices,1);
0260 
0261 f2 = zeros(Nvert,1);
0262 L  = zeros(Nvert,1);
0263 
0264 for index = 1:Nvert,
0265 
0266     % Find neighbouring vertex coordinates
0267     vi = find(FV.edge(index,:));  % the indices
0268 
0269     % Find distances between vertex and neighbours, using edge lengths.
0270     % The mean value is l in Smith (2002, eq. 4)
0271     edge_distance = FV.edge(index,vi);
0272     L(index) = mean(edge_distance);
0273 
0274     % Sn comes from update1
0275     sn = Sn(index,:);
0276     sn_mag = vector_magnitude(sn);
0277 
0278     % Calculate radius of local curvature, solve Smith (2002, eq. 4)
0279     radius_of_curvature = (L(index) ^ 2) / (2 * sn_mag);
0280 
0281     f2(index) = (1 + tanh( F * (1 / radius_of_curvature - E))) / 2;
0282 
0283 end
0284 
0285 return
0286 
0287 
0288 
0289 
0290 
0291 
0292 % --------------------------------------------------------
0293 %function [ f3, t1, Imin, Imax] = avw_bet_update3(FV,avw,bt,v,sn,t,t02,tm)
0294 function [ f3, t1, Imin, Imax] = avw_bet_update3(FV,avw,bt,sn,t,t02,tm)
0295 % component 3 of Smith (2002)
0296 
0297 
0298 %along a line pointing inward from the current vertex, min and max
0299 %intensities are found. d1 determines how far into the brain the minimum
0300 %intensity is searched for. d2 determines the same for the max intensity.
0301 
0302 d1 = 20; % 20 mm (Smith)
0303 d2 = d1/2;
0304 
0305 Nvert = size(FV.vertices,1);
0306 
0307 t1 = zeros(Nvert,1);
0308 f3 = zeros(Nvert,1);
0309 Imin = zeros(Nvert,1);
0310 Imax = zeros(Nvert,1);
0311 
0312 for index = 1:Nvert,
0313 
0314     x = FV.vertices(index,1);
0315     y = FV.vertices(index,2);
0316     z = FV.vertices(index,3);
0317     
0318     unit_normal = FV.unit_normals(index,:);
0319 
0320     % we need to find points along a line from the vertex point, parallel
0321     % with the vertex normal
0322 
0323     % vertex location (P) and unit vector of surface normal (n_hat), then
0324     % points (S) along line parallel to surface normal are obtained with vector
0325     % addition, such that S = P + (x * n_hat)
0326     d1_distances = [-d1:0.5:0 ]';
0327     d2_distances = [  0:0.5:d2]';
0328 
0329     XYZ = repmat([x y z],length(d1_distances),1);
0330     S1 = XYZ + (d1_distances * unit_normal);
0331     XYZ = repmat([x y z],length(d2_distances),1);
0332     S2 = XYZ + (d2_distances * unit_normal);
0333 
0334     X1 = S1(:,1);
0335     Y1 = S1(:,2);
0336     Z1 = S1(:,3);
0337 
0338     X2 = S2(:,1);
0339     Y2 = S2(:,2);
0340     Z2 = S2(:,3);
0341 
0342     % interpolate volume values at these points
0343     % ( not sure why have to swap XI,YI here )
0344     d1_intensity = interp3(avw.img,Y1,X1,Z1,'*nearest');
0345     d2_intensity = interp3(avw.img,Y2,X2,Z2,'*nearest');
0346         
0347     Imin(index) = max( [ t02; min([tm; d1_intensity]) ] );
0348     Imax(index) = min( [  tm; max([ t; d2_intensity]) ] );
0349 
0350     % for certain image intenisty distributions it can be varied (in the range
0351     % 0-1 to give optimal results, but the necessity for this is rare.
0352 
0353     t1(index) = (Imax(index) - t02) * bt + t02;
0354 
0355     f3(index) = ( 2 * (Imin(index) - t1(index)) ) / ( Imax(index) - t02 );
0356 
0357 end
0358 
0359 return
0360 
0361 
0362 
0363 
0364 
0365 
0366 
0367 
0368 
0369 
0370 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0371 % Radial update component 2, alternative to surface normal
0372 %
0373 %------------------------------------------------------
0374 function [FV,f2,sn,sn_unit,L] = avw_bet_update2_radial(FV,origin),
0375 
0376 % This function adapts Smith, S. (2002), Fast robust automated brain
0377 % extraction.  Human Brain Mapping, 17: 143-155.  This function
0378 % corresponds to update component 2: surface smoothness control.
0379 
0380 xo = origin(1); yo = origin(2); zo = origin(3);
0381 
0382 Nvert = size(FV.vertices,1);
0383 
0384 f2 = zeros(Nvert,1);
0385 sn = zeros(Nvert,3);
0386 
0387 for index = 1:Nvert,
0388 
0389     v = FV.vertices(index,:);
0390     x = FV.vertices(index,1);
0391     y = FV.vertices(index,2);
0392     z = FV.vertices(index,3);
0393 
0394     % Find radial distance of vertex from origin
0395     r = sqrt( (x-xo)2 + (y-yo)2 + (z-zo)2 );
0396 
0397     % Calculate unit vector
0398     v_unit_vector = ( v - origin ) / r;
0399 
0400     % Find directional cosines for line from center to vertex
0401     l = (x-xo)/r; % cos alpha
0402     m = (y-yo)/r; % cos beta
0403     n = (z-zo)/r; % cos gamma
0404 
0405     % Find neighbouring vertex coordinates
0406     vi = find(FV.edge(index,:));  % the indices
0407     neighbour_vertices = FV.vertices(vi,:);
0408     X = neighbour_vertices(:,1);
0409     Y = neighbour_vertices(:,2);
0410     Z = neighbour_vertices(:,3);
0411 
0412     % Find neighbour radial distances
0413     r_neighbours = sqrt( (X-xo).^2 + (Y-yo).^2 + (Z-zo).^2 );
0414     r_neighbours_mean = mean(r_neighbours);
0415     L(index) = r_neighbours_mean;
0416 
0417     % Find difference in radial distance between the vertex of interest and its
0418     % neighbours; this value approximates the magnitude of sn in
0419     % Smith (2002, eq. 1 to 4)
0420     r_diff = r - r_neighbours_mean;
0421 
0422     % Find the vector sn, in the direction of the vertex of interest, given the
0423     % difference in radial distance between vertex and mean of neighbours
0424 
0425     sn(index,:) = r_diff * v_unit_vector;
0426 
0427     snx = sn(index,1);
0428     sny = sn(index,2);
0429     snz = sn(index,3);
0430 
0431     sn_unit(index,:) = sn(index,:) ./ sqrt( (snx - v(1)).^2 + (sny - v(2)).^2 + (snz - v(3)).^2 );
0432 
0433     % Find distances between vertex and neighbours, using edge lengths.
0434     % The mean value is l in Smith (2002, eq. 4)
0435     edge_distance = FV.edge(index,vi);
0436     edge_distance_mean = mean(edge_distance);
0437 
0438     % Calculate radius of local curvature, solve Smith (2002, eq. 4)
0439     if r_diff,
0440         radius_of_curvature = (edge_distance_mean ^ 2) / (2 * r_diff);
0441     else
0442         radius_of_curvature = 10000;
0443     end
0444 
0445     % Define limits for radius of curvature
0446     radius_min =  3.33; % mm
0447     radius_max = 10.00; % mm
0448 
0449     % Sigmoid function parameters,
0450     % "where E and F control the scale and offset of the sigmoid"
0451     E = mean([(1 / radius_min),  (1 / radius_max)]);
0452     F = 6 * ( (1 / radius_min) - (1 / radius_max) );
0453 
0454     f2(index) = (1 + tanh( F * (1 / radius_of_curvature - E))) / 2;
0455 
0456 end
0457 
0458 return