Home > scripts > life_demo.m

life_demo

PURPOSE ^

% Example of initialization and fitting of the LiFE model

SYNOPSIS ^

function [fh, fe] = life_demo()

DESCRIPTION ^

% Example of initialization and fitting of the LiFE model
 This demo function illustrates how to:
  - A - Set up a LiFE structure, identified as 'fe' (fascicle evaluation) in
  the code below. This model contains a prediction of the diffusion
  measurements in each white-matter voxel given by the fascicles contained
  in a tractogrpahy solution, the connectome. Each fascicles makes a
  prediction about the direction of diffusion in the set of voxels where
  it travels through. The prediction is generated given the fascicle
  orientation and position isndie the voxel. Predictions from multiple
  fascicles in in each voxels are combined to generate a global connectome
  prediciton for the diffusion signal in large sets of white matter
  voxels.
  - B - Fit the LiFE model to compute the weights associated to each fascicle
  in the connectome. Fascicles in the conenctome contribute differently to
  predicting the diffusion signal in each voxel. First of all, fascicles
  make predictions about the diffusion only in voxels where they travel.
  Secondly, some fascicles have paths that produce better diffusion
  predictions than others. We use a least-square method to find the
  contribution of each fascicle to the diffusion signal in the white
  matter voxels where the fascicles travels. A single weight is assigned
  to each fascicle representing the global contribution of the fasicle to
  the signal of all the voxels along its path - we call this
  fascicle-global. Because multiple fascicles exist in several voxels the
  set of fascicles weights and fascicles predicitons represents the
  connectome-global prediction of the diffusion signal within the entire
  set of white matter voxels. Estimating the fascicle weights allows for
  evaluating the quality of the tractography solution. Eliminating
  fascicles that do not contribute to predicting the diffusion signal
  (they have assigned a zreo-weight). Finaly, the root-mean-squared error
  (RMSE) of the model to the diffusion data - the model prediction error -
  is used to evaluate the model prediction quality, compare different
  tractography models and to perform statistical inference on the on
  properties of the connectomes.
  - C - Compare two different connectome models. This demo will show how to
  compare two different conenctome models by using the diffusion
  prediction error (the Root-Mean-Squared Error, RMSE). We report the
  example of two conenctomes one generated using Constrained-spherical
  deconvolution (CSD) and probabilistic tractography the other using a
  tensor model and deterministic tractography
  - Note - The example connectomes used for this demo comprise a portion
  of the right occiptial lobe of an individual human brain. LiFE utilizes
  large-scale methods to solve the foward model. The software allows for
  solving connectomes spanning the entire white-matter of idnvidual
  brains. The size of the connectome on the test data set is small enought
  to allow for testing the code within a few minutes requiring only about
  10GB of computer RAM and standard hardaware. This code has been tested
  with: 

   - Ubuntu 12.04.4 LTS (Precise) 
   - 2.6Ghz i7 Processor and 24GB of RAM.
   - MatLab Version: 8.0.0.783 (R2012b)

   - Mac OSX 10.9
   - 1.8 Ghz i5 Processor and 8GB of RAM
   - MatLab Version: 8.0.0.783 (R2012b)

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Wed 16-Jul-2014 19:56:13 by m2html © 2005