feConnectomeCull

 Fits a connectome over and over by removing the fibers that make no
 contribution to the diffusion signal at each iteration.

   fe = feConnectomeCull(fe)

 Copyright (2013-2014), Franco Pestilli, Stanford University, pestillifranco@gmail.com.

0001 function [fe, cull] = feConnectomeCull(fe,display)
0002 % Fits a connectome over and over by removing the fibers that make no
0003 % contribution to the diffusion signal at each iteration.
0004 %
0005 %   fe = feConnectomeCull(fe)
0006 %
0007 % Copyright (2013-2014), Franco Pestilli, Stanford University, pestillifranco@gmail.com.
0008 
0009 if notDefined('display'); display.cull = false;end
0010 
0011 % Initialize a structure containing the information abotu the culling
0012 cull.name = feGet(fe,'name');
0013 cull.iter = 1;
0014 cull.minwval = eps; % (here we could have used 0 as thresholds, we use eps instead
0015                     % in hopes to seepd things up)
0016 
0017 % Extract the connectome model and the signa for fitting.
0018 M    = feGet(fe,'mfiber');
0019 dSig = feGet(fe,'dsigdemeaned');
0020 fe   = feSet(fe,'fit',feFitModel(M,dSig,'bbnnls'));
0021 
0022 % Count the number of fbers that make no contribution to the diffusion
0023 % prediction (here we could have used 0 as thresholds, we use eps instead
0024 % in hopes to seepd things up)
0025 sw   = sum(feGet(fe,'fiber weights') <= cull.minwval);
0026 
0027 % Update the culling structure given the results of the firts fit.
0028 cull.num2delete(cull.iter) = sw;
0029 cull.numtotal(cull.iter)   = size(M,2);
0030 cull.rmse(cull.iter)       = mean(feGet(fe,'vox rmse'));
0031 cull.rmse(cull.iter)       = mean(feGetRep(fe,'vox rmse'));
0032 cull.rrmse(cull.iter)      = mean(feGetRep(fe,'vox rmse ratio'));
0033 
0034 % Cull the connectome. We stop when no fibers have a contribution to the
0035 % diffusion signal less than eps.
0036 while (sw ~= 0)
0037     fprintf('Number of total %i and zero-weight %i fibers | culling...\n', ...
0038         cull.numtotal(cull.iter),cull.num2delete(cull.iter))
0039     M    = feGet(fe,'mfiber');
0040     dSig = feGet(fe,'dsigdemeaned');
0041     fe   = feSet(fe,'fit',feFitModel(M,dSig,'bbnnls'));
0042     sw   = sum(feGet(fe,'fiber weights') <= cull.minwval);
0043     fe   = feConnectomeReduceFibers(fe, find((feGet(fe,'fiber weights') > cull.minwval)));
0044     cull.iter = cull.iter + 1;
0045     cull.num2delete(cull.iter) = sw;
0046     cull.numtotal(cull.iter)   = size(M,2);
0047     cull.rmse(cull.iter)       = mean(feGet(fe,'vox rmse'));
0048     cull.rmse(cull.iter)       = mean(feGetRep(fe,'vox rmse'));
0049     cull.rrmse(cull.iter)      = mean(feGetRep(fe,'vox rmse ratio'));
0050 end
0051 
0052 % If display cull process
0053 if display.cull
0054 keyboard    
0055 end
0056 
0057 end