s_pestilli_etal_figures_4_5

% Example of initialization and fitting of the LiFE model

 This function illustrates how:
  - Set up an fe (fascicle evaluation) structure.
  - Combines them with fibers tracted between the two ROIs 
  - Generates an optimized connectome from a candidate connectome using LIFE 
  - Performs a virtual lesion

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

0001 function [fh, fe] = s_pestilli_etal_Figures_4_5()
0002 %% Example of initialization and fitting of the LiFE model
0003 %
0004 % This function illustrates how:
0005 %  - Set up an fe (fascicle evaluation) structure.
0006 %  - Combines them with fibers tracted between the two ROIs
0007 %  - Generates an optimized connectome from a candidate connectome using LIFE
0008 %  - Performs a virtual lesion
0009 %
0010 % Copyright (2013-2014), Franco Pestilli, Stanford University, pestillifranco@gmail.com.
0011 
0012 % Intialize a local matlab cluster if the parallel toolbox is available.
0013 feOpenLocalCluster;
0014 
0015 %% Build the file names for the diffusion data, the anatomical MR.
0016 dwiFile       = fullfile(lifeDemoDataPath('diffusion'),'pestilli_etal_life_demo_scan1_subject1_b2000_150dirs_stanford.nii.gz');
0017 dwiFileRepeat = fullfile(lifeDemoDataPath('diffusion'),'pestilli_etal_life_demo_scan2_subject1_b2000_150dirs_stanford.nii.gz');
0018 t1File        = fullfile(lifeDemoDataPath('anatomy'),  'pestilli_etal_life_demo_anatomy_t1w_stanford.nii.gz');
0019 
0020 %% (1) Probabilistic CSD-based connectome:
0021 % We will analyze first the CSD-based probabilistic tractography
0022 % connectome.
0023 fgFileName    = fullfile(lifeDemoDataPath('tractography'), ...
0024                 'pestilli_et_al_life_demo_mrtrix_csd_lmax10_probabilistic.mat');
0025 
0026 % The final connectome and data astructure will be saved with this name:
0027 feFileName    = 'life_build_model_demo_CSD_PROB';
0028 
0029 
0030 % Initialize the Connectome.
0031 fe = feConnectomeInit(dwiFile,fgFileName,feFileName,[],dwiFileRepeat,t1File);
0032 
0033 % Fit the model.
0034 fe = feSet(fe,'fit',feFitModel(feGet(fe,'mfiber'),feGet(fe,'dsigdemeaned'),'bbnnls'));
0035 
0036 % Plot a histogram of the RMSE.
0037 [fh(1), ~, ~] = plotHistRMSE(fe,'Probabilistic');
0038 
0039 % Plot a histogram of the RMSE ratio.
0040 % The Rrmse is the ratio between data test-retest reliability and
0041 % model error (the quality of the model fit).
0042 [fh(2), ~] = plotHistRrmse(fe,'Probabilistic');
0043 
0044 % Plot a histogram of the fitted fascicle weights.
0045 [fh(3), ~] = plotHistWeigths(fe,'Probabilistic');
0046 
0047 %% (2) Deterministic tensor-based connectome:
0048 % We will analyze first the CSD-based probabilistic tractography
0049 % connectome.
0050 fgFileName    = fullfile(lifeDemoDataPath('tractography'), ...
0051                 'pestilli_et_al_life_demo_mrtrix_tensor_deterministic.mat');
0052 
0053 % The final connectome and data astructure will be saved with this name:
0054 feFileName    = 'life_build_model_demo_TENSOR_DET';
0055 
0056 % Initialize the Connectome.
0057 fe = feConnectomeInit(dwiFile,fgFileName,feFileName,[],dwiFileRepeat,t1File);
0058 
0059 % Fit the model.
0060 fe = feSet(fe,'fit',feFitModel(feGet(fe,'mfiber'),feGet(fe,'dsigdemeaned'),'bbnnls'));
0061 
0062 % Plot a histogram of the RMSE.
0063 [fh(1), ~, ~] = plotHistRMSE(fe,'Deterministic');
0064 
0065 % Plot a histogram of the RMSE ratio.
0066 % The Rrmse is the ratio between data test-retest reliability and
0067 % model error (the quality of the model fit).
0068 [fh(2), ~] = plotHistRrmse(fe,'Deterministic');
0069 
0070 % Plot a histogram of the fitted fascicle weights.
0071 [fh(3), ~] = plotHistWeigths(fe,'Deterministic');
0072 
0073 keyboard
0074 end
0075 
0076 % ---------- Local  Functions ----------- %
0077 function [fh, rmse, rmsexv] = plotHistRMSE(fe,tractograpy)
0078 % Make a plot of the RMSE:
0079 
0080 % Extract the RMSE of the model on the fitted data set
0081 rmse   = feGet(fe,'vox rmse');
0082 
0083 % Extract the RMSE of the model on the second data set (cross-validated
0084 % error)
0085 rmsexv = feGetRep(fe,'vox rmse');
0086 
0087 figName = sprintf('%s - RMSE',tractograpy);
0088 fh = mrvNewGraphWin(figName);
0089 [y,x] = hist(rmse,50);
0090 plot(x,y,'k-');
0091 hold on
0092 [y,x] = hist(rmsexv,50);
0093 plot(x,y,'r-');
0094 set(gca,'tickdir','out','fontsize',16,'box','off');
0095 title('Root-mean squared error distribution across voxels','fontsize',16);
0096 ylabel('number of voxels','fontsize',16);
0097 xlabel('rmse (scanner units)','fontsize',16);
0098 legend({'RMSE fitted data set','RMSE cross-validated'},'fontsize',16);
0099 end
0100 
0101 function [fh, R] = plotHistRrmse(fe,tractograpy)
0102 % Make a plot of the RMSE Ratio:
0103 R       = feGetRep(fe,'voxrmseratio');
0104 figName = sprintf('%s - RMSE RATIO',tractograpy);
0105 fh      = mrvNewGraphWin(figName);
0106 [y,x]   = hist(R,linspace(.5,4,50));
0107 plot(x,y,'k-','linewidth',2);
0108 hold on
0109 plot([median(R) median(R)],[0 1200],'r-','linewidth',2);
0110 plot([1 1],[0 1200],'k-');
0111 set(gca,'tickdir','out','fontsize',16,'box','off');
0112 title('Root-mean squared error ratio','fontsize',16);
0113 ylabel('number of voxels','fontsize',16);
0114 xlabel('R_{rmse}','fontsize',16);
0115 legend({sprintf('Distribution of R_{rmse}'),sprintf('Median R_{rmse}')});
0116 end
0117 
0118 function [fh, w] = plotHistWeigths(fe,tractograpy)
0119 % Make a plot of the weights:
0120 w       = feGet(fe,'fiber weights');
0121 figName = sprintf('%s - Distribution of fascicle weights',tractograpy);
0122 fh      = mrvNewGraphWin(figName);
0123 [y,x]   = hist(w( w > 0 ),logspace(-5,-.3,40));
0124 semilogx(x,y,'k-','linewidth',2)
0125 set(gca,'tickdir','out','fontsize',16,'box','off')
0126 title( ...
0127     sprintf('Number of fascicles candidate connectome: %2.0f\nNumber of fascicles in optimized connetome: %2.0f' ...
0128     ,length(w),sum(w > 0)),'fontsize',16)
0129 ylabel('Number of fascicles','fontsize',16)
0130 xlabel('Fascicle weight','fontsize',16)
0131 end
0132