Description of life

0001 function [fh, fe] = life_demo()
0002 %% Example of initialization and fitting of the LiFE model
0003 % This demo function illustrates how to:
0004 %  - A - Set up a LiFE structure, identified as 'fe' (fascicle evaluation) in
0005 %  the code below. This model contains a prediction of the diffusion
0006 %  measurements in each white-matter voxel given by the fascicles contained
0007 %  in a tractogrpahy solution, the connectome. Each fascicles makes a
0008 %  prediction about the direction of diffusion in the set of voxels where
0009 %  it travels through. The prediction is generated given the fascicle
0010 %  orientation and position isndie the voxel. Predictions from multiple
0011 %  fascicles in in each voxels are combined to generate a global connectome
0012 %  prediciton for the diffusion signal in large sets of white matter
0013 %  voxels.
0014 %  - B - Fit the LiFE model to compute the weights associated to each fascicle
0015 %  in the connectome. Fascicles in the conenctome contribute differently to
0016 %  predicting the diffusion signal in each voxel. First of all, fascicles
0017 %  make predictions about the diffusion only in voxels where they travel.
0018 %  Secondly, some fascicles have paths that produce better diffusion
0019 %  predictions than others. We use a least-square method to find the
0020 %  contribution of each fascicle to the diffusion signal in the white
0021 %  matter voxels where the fascicles travels. A single weight is assigned
0022 %  to each fascicle representing the global contribution of the fasicle to
0023 %  the signal of all the voxels along its path - we call this
0024 %  fascicle-global. Because multiple fascicles exist in several voxels the
0025 %  set of fascicles weights and fascicles predicitons represents the
0026 %  connectome-global prediction of the diffusion signal within the entire
0027 %  set of white matter voxels. Estimating the fascicle weights allows for
0028 %  evaluating the quality of the tractography solution. Eliminating
0029 %  fascicles that do not contribute to predicting the diffusion signal
0030 %  (they have assigned a zreo-weight). Finaly, the root-mean-squared error
0031 %  (RMSE) of the model to the diffusion data - the model prediction error -
0032 %  is used to evaluate the model prediction quality, compare different
0033 %  tractography models and to perform statistical inference on the on
0034 %  properties of the connectomes.
0035 %  - C - Compare two different connectome models. This demo will show how to
0036 %  compare two different conenctome models by using the diffusion
0037 %  prediction error (the Root-Mean-Squared Error, RMSE). We report the
0038 %  example of two conenctomes one generated using Constrained-spherical
0039 %  deconvolution (CSD) and probabilistic tractography the other using a
0040 %  tensor model and deterministic tractography
0041 %  - Note - The example connectomes used for this demo comprise a portion
0042 %  of the right occiptial lobe of an individual human brain. LiFE utilizes
0043 %  large-scale methods to solve the foward model. The software allows for
0044 %  solving connectomes spanning the entire white-matter of idnvidual
0045 %  brains. The size of the connectome on the test data set is small enought
0046 %  to allow for testing the code within a few minutes requiring only about
0047 %  10GB of computer RAM and standard hardaware. This code has been tested
0048 %  with:
0049 %
0050 %   - Ubuntu 12.04.4 LTS (Precise)
0051 %   - 2.6Ghz i7 Processor and 24GB of RAM.
0052 %   - MatLab Version: 8.0.0.783 (R2012b)
0053 %
0054 %   - Mac OSX 10.9
0055 %   - 1.8 Ghz i5 Processor and 8GB of RAM
0056 %   - MatLab Version: 8.0.0.783 (R2012b)
0057 %
0058 % Copyright (2013-2014), Franco Pestilli, Stanford University, pestillifranco@gmail.com.
0059 
0060 % Intialize a local matlab cluster if the parallel toolbox is available.
0061 % This helps speeding up computations espacially for large conenctomes.
0062 feOpenLocalCluster;
0063 
0064 %% Build the file names for the diffusion data, the anatomical MRI.
0065 dwiFile       = fullfile(lifeDemoDataPath('diffusion'),'life_demo_scan1_subject1_b2000_150dirs_stanford.nii.gz');
0066 dwiFileRepeat = fullfile(lifeDemoDataPath('diffusion'),'life_demo_scan2_subject1_b2000_150dirs_stanford.nii.gz');
0067 t1File        = fullfile(lifeDemoDataPath('anatomy'),  'life_demo_anatomy_t1w_stanford.nii.gz');
0068 
0069 %% (1) Evaluate the Probabilistic CSD-based connectome.
0070 % We will analyze first the CSD-based probabilistic tractography
0071 % connectome.
0072 prob.tractography = 'Probabilistic';
0073 fgFileName    = fullfile(lifeDemoDataPath('tractography'), ...
0074                 'life_demo_mrtrix_csd_lmax10_probabilistic.mat');
0075 
0076 % The final connectome and data astructure will be saved with this name:
0077 feFileName    = 'life_build_model_demo_CSD_PROB';
0078 
0079 %% (1.1) Initialize the LiFE model structure, 'fe' in the code below.
0080 % This structure contains the forward model of diffusion based on the
0081 % tractography solution. It also contains all the information necessary to
0082 % compute model accuracry, and perform statistical tests. You can type
0083 % help('feBuildModel') in the MatLab prompt for more information.
0084 fe = feConnectomeInit(dwiFile,fgFileName,feFileName,[],dwiFileRepeat,t1File);
0085 
0086 %% (1.2) Fit the model.
0087 % Hereafter we fit the forward model of tracrography using a least-squared
0088 % method. The information generated by fitting the model (fiber weights
0089 % etc) is then installed in the LiFE structure.
0090 fe = feSet(fe,'fit',feFitModel(feGet(fe,'mfiber'),feGet(fe,'dsigdemeaned'),'bbnnls'));
0091 
0092 %% (1.3) Extract the RMSE of the model on the fitted data set.
0093 % We now use the LiFE structure and the fit to compute the error in each
0094 % white-matter voxel spanned by the tractography model.
0095 prob.rmse   = feGet(fe,'vox rmse');
0096 
0097 %% (1.4) Extract the RMSE of the model on the second data set.
0098 % Here we show how to compute the cross-valdiated RMSE of the tractography
0099 % model in each white-matter voxel. We store this information for later use
0100 % and to save computer memory.
0101 prob.rmsexv = feGetRep(fe,'vox rmse');
0102 
0103 %% (1.5) Extract the Rrmse.
0104 % We show how to extract the ratio between the model prediction error
0105 % (RMSE) and the test-retest reliability of the data.
0106 prob.rrmse  = feGetRep(fe,'vox rmse ratio');
0107 
0108 %% (1.6) Extract the fitted weights for the fascicles.
0109 % The following line shows how to extract the weight assigned to each
0110 % fascicle in the connectome.
0111 prob.w      = feGet(fe,'fiber weights');
0112 
0113 %% (1.7) Plot a histogram of the RMSE.
0114 % We plot the histogram of  RMSE across white-mater voxels.
0115 [fh(1), ~, ~] = plotHistRMSE(prob);
0116 
0117 %% (1.8) Plot a histogram of the RMSE ratio.
0118 % As a reminder the Rrmse is the ratio between data test-retest reliability
0119 % and model error (the quality of the model fit).
0120 [fh(2), ~] = plotHistRrmse(prob);
0121 
0122 %% (1.9) Plot a histogram of the fitted fascicle weights.
0123 [fh(3), ~] = plotHistWeights(prob);
0124 fe = feConnectomeInit(dwiFile,fgFileName,feFileName,[],dwiFileRepeat,t1File);
0125 
0126 %% Extract the coordinates of the white-matter voxels
0127 % We will use this later to compare probabilistic and deterministic models.
0128 p.coords = feGet(fe,'roi coords');
0129 clear fe
0130 
0131 %% (2) Evaluate the Deterministic tensor-based connectome.
0132 % We will now analyze the tensor-based Deterministic tractography
0133 % connectome.
0134 det.tractography = 'Deterministic';
0135 fgFileName    = fullfile(lifeDemoDataPath('tractography'), ...
0136                 'life_demo_mrtrix_tensor_deterministic.mat');
0137 
0138 % The final connectome and data astructure will be saved with this name:
0139 feFileName    = 'life_build_model_demo_TENSOR_DET';
0140 
0141 %% (2.1) Initialize the LiFE model structure, 'fe' in the code below.
0142 % This structure contains the forward model of diffusion based on the
0143 % tractography solution. It also contains all the information necessary to
0144 % compute model accuracry, and perform statistical tests. You can type
0145 % help('feBuildModel') in the MatLab prompt for more information.
0146 fe = feConnectomeInit(dwiFile,fgFileName,feFileName,[],dwiFileRepeat,t1File);
0147 
0148 %% (2.2) Fit the model.
0149 % Hereafter we fit the forward model of tracrography using a least-squared
0150 % method. The information generated by fitting the model (fiber weights
0151 % etc) is then installed in the LiFE structure.
0152 fe = feSet(fe,'fit',feFitModel(feGet(fe,'mfiber'),feGet(fe,'dsigdemeaned'),'bbnnls'));
0153 
0154 %% (2.3) Extract the RMSE of the model on the fitted data set.
0155 % We now use the LiFE structure and the fit to compute the error in each
0156 % white-matter voxel spanned by the tractography model.
0157 det.rmse   = feGet(fe,'vox rmse');
0158 
0159 %% (2.4) Extract the RMSE of the model on the second data set.
0160 % Here we show how to compute the cross-valdiated RMSE of the tractography
0161 % model in each white-matter voxel. We store this information for later use
0162 % and to save computer memory.
0163 det.rmsexv = feGetRep(fe,'vox rmse');
0164 
0165 %% (2.5) Extract the Rrmse.
0166 % We show how to extract the ratio between the model prediction error
0167 % (RMSE) and the test-retest reliability of the data.
0168 det.rrmse  = feGetRep(fe,'vox rmse ratio');
0169 
0170 %% (2.6) Extract the fitted weights for the fascicles.
0171 % The following line shows how to extract the weight assigned to each
0172 % fascicle in the connectome.
0173 det.w      = feGet(fe,'fiber weights');
0174 
0175 %% (2.7) Plot a histogram of the RMSE.
0176 % We plot the histogram of  RMSE across white-mater voxels.
0177 [fh(1), ~, ~] = plotHistRMSE(det);
0178 
0179 %% (2.8) Plot a histogram of the RMSE ratio.
0180 % As a reminder the Rrmse is the ratio between data test-retest reliability
0181 % and model error (the quality of the model fit).
0182 [fh(2), ~] = plotHistRrmse(det);
0183 
0184 %% (2.9) Plot a histogram of the fitted fascicle weights.
0185 [fh(3), ~] = plotHistWeights(det);
0186 
0187 %% Extract the coordinates of the white-matter voxels.
0188 % We will use this later to compare probabilistic and deterministic models.
0189 d.coords = feGet( fe, 'roi coords');
0190 clear fe
0191 
0192 %% (3) Compare the quality of fit of Probabilistic and Deterministic connectomes.
0193 %% (3.1) Find the common coordinates between the two connectomes.
0194 %
0195 % The two tractography method might have passed through slightly different
0196 % white-matter voxels. Here we find the voxels where both models passed. We
0197 % will compare the error only in these common voxels. There are more
0198 % coordinates in the Prob connectome, because the tracking fills up more
0199 % White-matter.
0200 %
0201 % So, hereafter:
0202 % - First we find the indices in the probabilistic connectome of the
0203 % coordinate in the deterministic connectome. But there are some of the
0204 % coordinates in the Deterministic conectome that are NOT in the
0205 % Probabilistic connectome.
0206 %
0207 % - Second we find the indices in the Deterministic connectome of the
0208 % subset of coordinates in the Probabilistic connectome found in the
0209 % previous step.
0210 %
0211 % - Third we find the common voxels. These allow us to find the rmse for
0212 % the same voxels.
0213 fprintf('Finding common brain coordinates between P and D connectomes...\n')
0214 prob.coordsIdx = ismember(p.coords,d.coords,'rows');
0215 prob.coords    = p.coords(prob.coordsIdx,:);
0216 det.coordsIdx  = ismember(d.coords,prob.coords,'rows');
0217 det.coords     = d.coords(det.coordsIdx,:);
0218 prob.rmse      = prob.rmse( prob.coordsIdx);
0219 det.rmse       = det.rmse( det.coordsIdx);
0220 clear p d
0221 
0222 %% (3.2) Make a scatter plot of the RMSE of the two tractography models
0223 fh(4) = scatterPlotRMSE(det,prob);
0224 
0225 %% (3.3) Compute the strength-of-evidence (S) and the Earth Movers Distance.
0226 % Compare the RMSE of the two models using the Stregth-of-evidence and the
0227 % Earth Movers Distance.
0228 se = feComputeEvidence(prob.rmse,det.rmse);
0229 
0230 %% (3.4) Strength of evidence in favor of Probabilistic tractography.
0231 % Plot the distributions of resampled mean RMSE
0232 % used to compute the strength of evidence (S).
0233 fh(5) = distributionPlotStrengthOfEvidence(se);
0234 
0235 %% (3.5) RMSE distributions for Probabilistic and Deterministic tractography.
0236 % Compare the distributions using the Earth Movers Distance.
0237 % Plot the distributions of RMSE for the two models and report the Earth
0238 % Movers Distance between the distributions.
0239 fh(6) = distributionPlotEarthMoversDistance(se);
0240 
0241 end
0242 
0243 % ---------- Local Plot Functions ----------- %
0244 function [fh, rmse, rmsexv] = plotHistRMSE(info)
0245 % Make a plot of the RMSE:
0246 rmse   = info.rmse;
0247 rmsexv = info.rmsexv;
0248 
0249 figName = sprintf('%s - RMSE',info.tractography);
0250 fh = mrvNewGraphWin(figName);
0251 [y,x] = hist(rmse,50);
0252 plot(x,y,'k-');
0253 hold on
0254 [y,x] = hist(rmsexv,50);
0255 plot(x,y,'r-');
0256 set(gca,'tickdir','out','fontsize',16,'box','off');
0257 title('Root-mean squared error distribution across voxels','fontsize',16);
0258 ylabel('number of voxels','fontsize',16);
0259 xlabel('rmse (scanner units)','fontsize',16);
0260 legend({'RMSE fitted data set','RMSE cross-validated'},'fontsize',16);
0261 end
0262 
0263 function [fh, R] = plotHistRrmse(info)
0264 % Make a plot of the RMSE Ratio:
0265 
0266 R       = info.rrmse;
0267 figName = sprintf('%s - RMSE RATIO',info.tractography);
0268 fh      = mrvNewGraphWin(figName);
0269 [y,x]   = hist(R,linspace(.5,4,50));
0270 plot(x,y,'k-','linewidth',2);
0271 hold on
0272 plot([median(R) median(R)],[0 1200],'r-','linewidth',2);
0273 plot([1 1],[0 1200],'k-');
0274 set(gca,'tickdir','out','fontsize',16,'box','off');
0275 title('Root-mean squared error ratio','fontsize',16);
0276 ylabel('number of voxels','fontsize',16);
0277 xlabel('R_{rmse}','fontsize',16);
0278 legend({sprintf('Distribution of R_{rmse}'),sprintf('Median R_{rmse}')});
0279 end
0280 
0281 function [fh, w] = plotHistWeights(info)
0282 % Make a plot of the weights:
0283 
0284 w       = info.w;
0285 figName = sprintf('%s - Distribution of fascicle weights',info.tractography);
0286 fh      = mrvNewGraphWin(figName);
0287 [y,x]   = hist(w( w > 0 ),logspace(-5,-.3,40));
0288 semilogx(x,y,'k-','linewidth',2)
0289 set(gca,'tickdir','out','fontsize',16,'box','off')
0290 title( ...
0291     sprintf('Number of fascicles candidate connectome: %2.0f\nNumber of fascicles in optimized connetome: %2.0f' ...
0292     ,length(w),sum(w > 0)),'fontsize',16)
0293 ylabel('Number of fascicles','fontsize',16)
0294 xlabel('Fascicle weight','fontsize',16)
0295 end
0296 
0297 function fh = scatterPlotRMSE(det,prob)
0298 figNameRmse = sprintf('prob_vs_det_rmse_common_voxels_map');
0299 fh = mrvNewGraphWin(figNameRmse);
0300 [ymap,x]  = hist3([det.rmse;prob.rmse]',{[10:1:70], [10:1:70]});
0301 ymap = ymap./length(prob.rmse);
0302 sh   = imagesc(flipud(log10(ymap)));
0303 cm   = colormap(flipud(hot)); view(0,90);
0304 axis('square')      
0305 set(gca, ...
0306     'xlim',[1 length(x{1})],...
0307     'ylim',[1 length(x{1})], ...
0308     'ytick',[1 (length(x{1})/2) length(x{1})], ...
0309     'xtick',[1 (length(x{1})/2) length(x{1})], ...
0310     'yticklabel',[x{1}(end) x{1}(round(end/2)) x{1}(1)], ...
0311     'xticklabel',[x{1}(1)   x{1}(round(end/2)) x{1}(end)], ...
0312     'tickdir','out','ticklen',[.025 .05],'box','off', ...
0313     'fontsize',16,'visible','on')
0314 hold on
0315 plot3([1 length(x{1})],[length(x{1}) 1],[max(ymap(:)) max(ymap(:))],'k-','linewidth',1)
0316 ylabel('Deterministic_{rmse}','fontsize',16)
0317 xlabel('Probabilistic_{rmse}','fontsize',16)
0318 cb = colorbar;
0319 tck = get(cb,'ytick');
0320 set(cb,'yTick',[min(tck)  mean(tck) max(tck)], ...
0321     'yTickLabel',round(1000*10.^[min(tck),...
0322     mean(tck), ...
0323     max(tck)])/1000, ...
0324     'tickdir','out','ticklen',[.025 .05],'box','on', ...
0325     'fontsize',16,'visible','on')
0326 end
0327 
0328 function fh = distributionPlotStrengthOfEvidence(se)
0329 
0330 y_e        = se.s.unlesioned_e;
0331 ywo_e      = se.s.lesioned_e;
0332 dprime     = se.s.mean;
0333 std_dprime = se.s.std;
0334 xhis       = se.s.unlesioned.xbins;
0335 woxhis     = se.s.lesioned.xbins;
0336 
0337 histcolor{1} = [0 0 0];
0338 histcolor{2} = [.95 .6 .5];
0339 figName = sprintf('Strength_of_Evidence_test_PROB_vs_DET_model_rmse_mean_HIST');
0340 fh = mrvNewGraphWin(figName);
0341 patch([xhis,xhis],y_e(:),histcolor{1},'FaceColor',histcolor{1},'EdgeColor',histcolor{1});
0342 hold on
0343 patch([woxhis,woxhis],ywo_e(:),histcolor{2},'FaceColor',histcolor{2},'EdgeColor',histcolor{2}); 
0344 set(gca,'tickdir','out', ...
0345         'box','off', ...
0346         'ticklen',[.025 .05], ...
0347         'ylim',[0 .2], ... 
0348         'xlim',[min(xhis) max(woxhis)], ...
0349         'xtick',[min(xhis) round(mean([xhis, woxhis])) max(woxhis)], ...
0350         'ytick',[0 .1 .2], ...
0351         'fontsize',16)
0352 ylabel('Probability','fontsize',16)
0353 xlabel('rmse','fontsize',16')
0354 
0355 title(sprintf('Strength of evidence:\n mean %2.3f - std %2.3f',dprime,std_dprime), ...
0356     'FontSize',16)
0357 legend({'Probabilistic','Deterministic'})
0358 end
0359 
0360 function fh = distributionPlotEarthMoversDistance(se)
0361 
0362 prob = se.nolesion;
0363 det  = se.lesion;
0364 em   = se.em;
0365 
0366 histcolor{1} = [0 0 0];
0367 histcolor{2} = [.95 .6 .5];
0368 figName = sprintf('EMD_PROB_DET_model_rmse_mean_HIST');
0369 fh = mrvNewGraphWin(figName);
0370 plot(prob.xhist,prob.hist,'r-','color',histcolor{1},'linewidth',4);
0371 hold on
0372 plot(det.xhist,det.hist,'r-','color',histcolor{2},'linewidth',4); 
0373 set(gca,'tickdir','out', ...
0374         'box','off', ...
0375         'ticklen',[.025 .05], ...
0376         'ylim',[0 .12], ... 
0377         'xlim',[0 95], ...
0378         'xtick',[0 45 90], ...
0379         'ytick',[0 .06 .12], ...
0380         'fontsize',16)
0381 ylabel('Proportion white-matter volume','fontsize',16)
0382 xlabel('RMSE (raw MRI scanner units)','fontsize',16')
0383 title(sprintf('Earth Movers Distance: %2.3f (raw scanner units)',em.mean),'FontSize',16)
0384 legend({'Probabilistic','Deterministic'})
0385 end
0386
life_demo

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
