s_pestilli_etal_figure_7

% Comparison between two tractography model.

  This function illustrates how to:
  - Compute the Room-Mean-Square-Error (RMSE) from precomputed LiFE structures.
  - Compare the RMSE of two different tractography models. In this example
    we show how to compare a Probabilistic and a Deterministic
    tractogrpahy model.
  - We show how to use the LiFE software to compute the "Strength of evidence" and
    the "Earth Movers Distance" to assess the difference in RMSE between
    two tractography models.

 The example shows that for this brain and tractography model the
 Probabilitic tractorgrpahy model generates smaller error; it predicts the
 diffusion measuments better than the Deterministic.
 
 Results similar to the ones in reproduced in this exmple are reported in
 Fig. 7 of Pestilli et al. UNDER REVIEW.

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

0001 function s_pestilli_etal_figure_7()
0002 %% Comparison between two tractography model.
0003 %
0004 %  This function illustrates how to:
0005 %  - Compute the Room-Mean-Square-Error (RMSE) from precomputed LiFE structures.
0006 %  - Compare the RMSE of two different tractography models. In this example
0007 %    we show how to compare a Probabilistic and a Deterministic
0008 %    tractogrpahy model.
0009 %  - We show how to use the LiFE software to compute the "Strength of evidence" and
0010 %    the "Earth Movers Distance" to assess the difference in RMSE between
0011 %    two tractography models.
0012 %
0013 % The example shows that for this brain and tractography model the
0014 % Probabilitic tractorgrpahy model generates smaller error; it predicts the
0015 % diffusion measuments better than the Deterministic.
0016 %
0017 % Results similar to the ones in reproduced in this exmple are reported in
0018 % Fig. 7 of Pestilli et al. UNDER REVIEW.
0019 %
0020 % Copyright (2013-2014), Franco Pestilli, Stanford University, pestillifranco@gmail.com
0021 
0022 % Get the base directory for the data
0023 datapath = pestilliDataPath;
0024 
0025 %% Load two pre-computed connectomes generated using Probabilistic and Deterministic tractography.
0026 %
0027 feProbFileName = 'subject1_life_culled_2mm_150dir_b2000_probabilistic_lmax8_diffModAx100Rd0.mat';
0028 feDetFileName  = 'subject1_life_culled_2mm_150dir_b2000_tensor_diffModAx100Rd0.mat';
0029 
0030 fprintf('Loading precomputed LiFE models for probabilistic (P) and deterministic (D) connectomes...\n')
0031 p = load(fullfile(datapath,'life_structures',feProbFileName));
0032 d = load(fullfile(datapath,'life_structures',feDetFileName)); 
0033 
0034 %% Extract the RMSE of the LiFE model for each connectome in each white-matter voxel.
0035 % We compute the RMSE for each individual voxel within the White-matter.
0036 fprintf('Computing Root-Mean-Square-Error (RMSE) for each brain voxel and tractography model..\n')
0037 p.rmse   = feGetRep(p.fe,'vox rmse');
0038 d.rmse   = feGetRep(d.fe, 'vox rmse');
0039 
0040 %% Find the common coordinates between the two connectomes.
0041 % The two tractography method might have passed through slightly different
0042 % voxels. Here we find the voxels where both models passed. We will compare
0043 % the error only in these common voxels. There are more coordinates in the
0044 % Prob connectome, because the tracking fills up more White-matter.
0045 %
0046 % So, hereafter:
0047 %
0048 % - First we find the indices in the probabilistic connectome of the
0049 % coordinate in the deterministic connectome. But there are some of the
0050 % coordinates in the Deterministic conectome that are NOT in the
0051 % Probabilistic connectome.
0052 %
0053 % - Second we find the indices in the Deterministic connectome of the
0054 % subset of coordinates in the Probabilistic connectome found in the
0055 % previous step.
0056 %
0057 % - Third we find the common voxels. These allow us to find the rmse for
0058 % the same voxels.
0059 
0060 fprintf('Finding common brain coordinates between P and D connectomes...\n')
0061 p.coords = feGet(   p.fe,'roi coords');
0062 d.coords = feGet(   d.fe, 'roi coords');
0063 prob.coordsIdx = ismember(p.coords,d.coords,'rows');
0064 prob.coords   = p.coords(prob.coordsIdx,:);
0065 det.coordsIdx = ismember(d.coords,prob.coords,'rows');
0066 det.coords    = d.coords(det.coordsIdx,:);
0067 prob.rmse  = p.rmse( prob.coordsIdx);
0068 det.rmse   = d.rmse( det.coordsIdx);
0069 
0070 
0071 %% Compare the RMSE of the Probabilistic and Deterministic models using a scatter-density plot.
0072 scatterPlotRMSE(det,prob)
0073 
0074 %% Compare the RMSE of the two models using the Stregth-of-evidence and the Earth Movers Distance.
0075 se = feComputeEvidence(p.rmse,d.rmse);
0076 
0077 %% Show the strength of evidence in favor of Probabilistic versus Deterministic tractography.
0078 % Plot the distributions of resampled mean RMSE used to compute the strength of
0079 % evidence (S).
0080 distributionPlotStrengthOfEvidence(se.s.unlesioned_e,se.s.lesioned_e,se.s.mean,se.s.std,se.s.unlesioned.xbins, se.s.lesioned.xbins)
0081 
0082 %% Show the RMSE distributions for Probabilistic deterministic tractography.
0083 %% Compare the distributions using the Earth Movers Distance.
0084 % Plot the distributions of RMSE for the two models and report the Earth
0085 % Movers Distance between the distributions
0086 
0087 distributionPlotEarthMoversDistance(se.nolesion,se.lesion,se.em)
0088 
0089 end % End main function
0090 
0091 %% Helper functions used for plotting
0092 function scatterPlotRMSE(det,prob)
0093 figNameRmse = sprintf('prob_vs_det_rmse_common_voxels_map');
0094 mrvNewGraphWin(figNameRmse);
0095 [ymap,x]  = hist3([det.rmse;prob.rmse]',{[10:1:70], [10:1:70]});
0096 ymap = ymap./length(prob.rmse);
0097 sh   = imagesc(flipud(log10(ymap)));
0098 cm   = colormap(flipud(hot)); view(0,90);
0099 axis('square')      
0100 set(gca, ...
0101     'xlim',[1 length(x{1})],...
0102     'ylim',[1 length(x{1})], ...
0103     'ytick',[1 (length(x{1})/2) length(x{1})], ...
0104     'xtick',[1 (length(x{1})/2) length(x{1})], ...
0105     'yticklabel',[x{1}(end) x{1}(round(end/2)) x{1}(1)], ...
0106     'xticklabel',[x{1}(1)   x{1}(round(end/2)) x{1}(end)], ...
0107     'tickdir','out','ticklen',[.025 .05],'box','off', ...
0108     'fontsize',16,'visible','on')
0109 hold on
0110 plot3([1 length(x{1})],[length(x{1}) 1],[max(ymap(:)) max(ymap(:))],'k-','linewidth',1)
0111 ylabel('Deterministic_{rmse}','fontsize',16)
0112 xlabel('Probabilistic_{rmse}','fontsize',16)
0113 cb = colorbar;
0114 tck = get(cb,'ytick');
0115 set(cb,'yTick',[min(tck)  mean(tck) max(tck)], ...
0116     'yTickLabel',round(1000*10.^[min(tck),...
0117     mean(tck), ...
0118     max(tck)])/1000, ...
0119     'tickdir','out','ticklen',[.025 .05],'box','on', ...
0120     'fontsize',16,'visible','on')
0121 end
0122 
0123 function distributionPlotStrengthOfEvidence(y_e,ywo_e,dprime,std_dprime,xhis,woxhis)
0124 histcolor{1} = [0 0 0];
0125 histcolor{2} = [.95 .6 .5];
0126 figName = sprintf('Strength_of_Evidence_test_PROB_vs_DET_model_rmse_mean_HIST');
0127 mrvNewGraphWin(figName);
0128 patch([xhis,xhis],y_e(:),histcolor{1},'FaceColor',histcolor{1},'EdgeColor',histcolor{1});
0129 hold on
0130 patch([woxhis,woxhis],ywo_e(:),histcolor{2},'FaceColor',histcolor{2},'EdgeColor',histcolor{2}); 
0131 set(gca,'tickdir','out', ...
0132         'box','off', ...
0133         'ticklen',[.025 .05], ...
0134         'ylim',[0 .2], ... 
0135         'xlim',[28 34], ...
0136         'xtick',[28 30 32 34], ...
0137         'ytick',[0 .1 .2], ...
0138         'fontsize',16)
0139 ylabel('Probability','fontsize',16)
0140 xlabel('rmse','fontsize',16')
0141 
0142 title(sprintf('Strength of evidence:\n mean %2.3f - std %2.3f',dprime,std_dprime), ...
0143     'FontSize',16)
0144 legend({'Probabilistic','Deterministic'})
0145 end
0146 
0147 
0148 function distributionPlotEarthMoversDistance(prob,det,em)
0149 histcolor{1} = [0 0 0];
0150 histcolor{2} = [.95 .6 .5];
0151 figName = sprintf('EMD_PROB_DET_model_rmse_mean_HIST');
0152 mrvNewGraphWin(figName);
0153 plot(prob.xhist,prob.hist,'r-','color',histcolor{1},'linewidth',4);
0154 hold on
0155 plot(det.xhist,det.hist,'r-','color',histcolor{2},'linewidth',4); 
0156 set(gca,'tickdir','out', ...
0157         'box','off', ...
0158         'ticklen',[.025 .05], ...
0159         'ylim',[0 .12], ... 
0160         'xlim',[0 95], ...
0161         'xtick',[0 45 90], ...
0162         'ytick',[0 .06 .12], ...
0163         'fontsize',16)
0164 ylabel('Proportion white-matter volume','fontsize',16)
0165 xlabel('RMSE (raw MRI scanner units)','fontsize',16')
0166 title(sprintf('Earth Movers Distance: %2.3f (raw scanner units)',em.mean),'FontSize',16)
0167 legend({'Probabilistic','Deterministic'})
0168 end