feVirtualLesion

 Perform a virtual lesion.
 - Remove a set of fascicles from a whole-brain connectome
 - Find the fascicle's path-neighborhood (the set of fascicles that share the same voxels)
 - Reduce the connectome to only the voxels of the set of fascicels that
   we are lesioning. These are the voxels that contain the signal and the
   error useful for the lesion test.
 - Compute a series of statistics on the lesioned connectome.

 
 [se, fh,  feLesion,  feNoLesion, newFascicleIndices, indicesFibersKept, commonCoords] = ...
           feVirtualLesion(feNoLesion, fascicleIndices, refitConnectome, display)

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

0001 function [se, fig,  feLesion,  feNoLesion, newFascicleIndices, indicesFibersKept, commonCoords] = ...
0002           feVirtualLesion(feNoLesion, fascicleIndices, display, refitConnectome)
0003 % Perform a virtual lesion.
0004 % - Remove a set of fascicles from a whole-brain connectome
0005 % - Find the fascicle's path-neighborhood (the set of fascicles that share the same voxels)
0006 % - Reduce the connectome to only the voxels of the set of fascicels that
0007 %   we are lesioning. These are the voxels that contain the signal and the
0008 %   error useful for the lesion test.
0009 % - Compute a series of statistics on the lesioned connectome.
0010 %
0011 %
0012 % [se, fh,  feLesion,  feNoLesion, newFascicleIndices, indicesFibersKept, commonCoords] = ...
0013 %           feVirtualLesion(feNoLesion, fascicleIndices, refitConnectome, display)
0014 %
0015 %
0016 % Copyright (2013-2014), Franco Pestilli, Stanford University, pestillifranco@gmail.com.
0017 
0018 if notDefined('display'),   
0019     display.tract = 0;
0020     display.distributions = 0; 
0021     display.evidence = 0; 
0022     fig = []; 
0023 end
0024 if notDefined('refitConnectome'), refitConnectome = 0;end
0025 
0026 feOpenLocalCluster
0027 tic,fprintf('\n[%s] Performing the lesion (removing fascicles from the connectome)... ',mfilename)
0028 if ~islogical( fascicleIndices )
0029     nfibers = feGet(feNoLesion,'nfibers');
0030     fascicles2keep = false(nfibers,1);
0031     fascicles2keep(fascicleIndices) = true;
0032 else
0033     fascicles2keep = fascicleIndices;
0034 end
0035 feLesion = feConnectomeReduceFibers(feNoLesion, fascicles2keep );toc
0036 
0037 % Extract the fascicle out of the fiber group.
0038 fas = fgExtract(feGet(feNoLesion,'fibers img'), fascicleIndices, 'keep' );
0039 fasAcpc = fgExtract(feGet(feNoLesion,'fibers acpc'), fascicleIndices, 'keep' );
0040 
0041 % Get the cordinates of the fascicle we just deleted. These coordinates are
0042 % contained in the full connectome. We want to fin the indices in the M
0043 % matrix to the the voxels of the fascicle and keep only those voxels.
0044 tic, fprintf('\n[%s] Finding the voxels of fascicle... ',mfilename)
0045 fasCoords    = fefgGet(fas,'unique image coords');
0046 allCoords    = feGet(feNoLesion,   'roi coords');
0047 commonCoords = ismember(allCoords, fasCoords,  'rows'); toc
0048 
0049 % Now: commonCoords contains the indices of the voxels to keep from feLesion,
0050 % these voxels are part of the connectome feNoLesion. So now we want to delete all
0051 % the rest of the voxels from feLesion and keep only the voxels where
0052 % the fascicle in feFP goes through.
0053 %
0054 % At the same time we keep all the fibers in the connectome that still pass
0055 % throught the left voxels in the ROI. We will use this subset of voxels
0056 % and fibers to test the quality of fit of the connectome at the location
0057 % where the feFN fascicle went through.
0058 tic,fprintf('\n[%s] Creating a lesioned connectome... ',mfilename)
0059 feLesion = feConnectomeReduceVoxels(feLesion,find(commonCoords));
0060 toc
0061 
0062 tic,fprintf('\n[%s] Creating an unleasioned path-neighborhood connectome... ',mfilename)
0063 [feNoLesion, indicesFibersKept]    = feConnectomeReduceVoxels(feNoLesion, find(commonCoords));
0064 toc
0065 
0066 tic,fprintf('\n[%s] Finding the pathneighborhood voxels... ',mfilename)
0067 % Here we return the indices of the fascicle in the newly resized connectome.
0068 newFascicleIndices = ismember(find(indicesFibersKept),...
0069                               find(fascicleIndices),'rows');
0070 toc
0071 
0072 % Fit the connectomes unlesioned the leftover fibers.
0073 if refitConnectome
0074     % Refit and install the fit.
0075     tic, fprintf('\n[%s] Fitting connectome lesioned the fascicle... ',mfilename)
0076     feLesion = feSet(feLesion,'fit', ...
0077         feFitModel(feGet(feLesion,'Mfiber'),feGet(feLesion,'dsigdemeaned'), ...
0078         'bbnnls'));toc
0079     
0080     tic, fprintf('\n[%s] Fitting connectome unlesioned the fascicle... ',mfilename)
0081     feNoLesion    = feSet(feNoLesion,   'fit', ...
0082         feFitModel(feGet(feNoLesion,'Mfiber'),   feGet(feNoLesion,'dsigdemeaned'), ...
0083         'bbnnls'));toc
0084 end
0085 se = feComputeEvidence((feGetRep(feNoLesion,'vox  rmse')),(feGetRep(feLesion, 'vox  rmse')));
0086    
0087 fig(1).name = sprintf('rmse_distributions_%s',mfilename);
0088 fig(1).h   = nan;
0089 fig(1).type = 'eps';
0090 if display.distributions
0091     % Raw RMSE distirbutions
0092     fig(1).name = sprintf('rmse_distributions_%s',mfilename);
0093     fig(1).h   = figure('name',fig(1).name,'color','w');
0094     fig(1).type = 'eps';
0095     set(fig(1).h,'Units','normalized','Position',[0.007 0.55  0.28 0.36]);
0096     plot(se.lesion.xhist,se.lesion.hist,'-','color', [.95 .45 .1],'linewidth',2); hold on
0097     plot(se.nolesion.xhist,se.nolesion.hist,'-','linewidth',2, 'color', [.1 .45 .95])
0098     plot([se.nolesion.rmse.mean,se.nolesion.rmse.mean], [0,0.12],'-','color',[.1 .45 .95] ) 
0099     plot([se.lesion.rmse.mean,se.lesion.rmse.mean], [0,0.12], '-', 'color',[.95 .45 .1])
0100     title(sprintf('mean RMSE\nno-lesion %2.3f | lesion %2.2f',se.nolesion.rmse.mean,se.lesion.rmse.mean),'fontsize',16)
0101     ylabel('Probability', 'fontsize',14);xlabel('RMSE', 'fontsize',14)
0102     legend({'Lesion','No lesion'},'fontsize',14);
0103     set(gca,'box','off','xtick',[0 round(se.xrange(2)/2) se.xrange(2)],'ytick',[0 .06 .12],'xlim',[0 se.xrange(2)],'ylim',[0 .125], ...
0104         'tickdir', 'out', 'ticklength', [0.025 0])
0105 
0106     % Plot the null distribution and the empirical difference
0107     % We reorganize the variables names hee below just to keep the plotting
0108     % code more compact.
0109     ywo_e = se.s.lesioned_e;
0110     y_e   = se.s.unlesioned_e;
0111     woxhis = se.s.lesioned.xbins;
0112     xhis   = se.s.unlesioned.xbins;
0113     min_x = se.s.min_x;
0114     max_x = se.s.max_x; 
0115     fig(2).name = sprintf('virtual_lesion_test_mean_rmse_hist_%s_%s',mfilename,feLesion.name);
0116     fig(2).h   = figure('name',fig(2).name,'color','w');  
0117     fig(2).type = 'eps';
0118     set(fig(2).h,'Units','normalized','Position',[0.007 0.55  0.28 0.36]);
0119     patch([xhis,xhis],y_e(:),[.1 .45 .95],'FaceColor',[.1 .45 .95],'EdgeColor',[.1 .45 .95]);
0120     hold on
0121     patch([woxhis,woxhis],ywo_e(:),[.95 .45 .1],'FaceColor',[.95 .45 .1],'EdgeColor',[.95 .45 .1]);
0122     set(gca,'tickdir','out', ...
0123         'box','off', ...
0124         'ylim',[0 0.25], ...
0125         'xlim',[min_x,max_x], ...
0126         'ytick',[0 0.1 0.2], ...
0127         'xtick',round(linspace(min_x,max_x,4)), ...
0128         'fontsize',16)
0129     ylabel('Probability','fontsize',16)
0130     xlabel('rmse','fontsize',16')
0131     title(sprintf('Strength of connection evidence %2.3f',(se.s.mean)), 'FontSize',16)
0132 end
0133 
0134 fig(3).name = sprintf('tract_%s',mfilename);
0135 fig(3).h   = nan;
0136 fig(3).type = 'jpg';
0137 fig(4).name = sprintf('pathneighborhood_%s',mfilename);
0138 fig(4).h   = nan;
0139 fig(4).type = 'jpg';
0140 fig(5).name = sprintf('pathneighborhood_AND_tract_%s',mfilename);
0141 fig(5).h   = nan;
0142 fig(5).type = 'jpg';
0143 if display.tract
0144     % Now display the lesion performed.
0145     % (1) The fascicle to delete.
0146     fig(3).name = sprintf('tract_%s',mfilename);
0147     fig(3).h   = figure('name',fig(3).name,'color',[.1 .45 .95]); 
0148     fig(3).type = 'jpg';
0149     [fig(3).h, fig(3).light] =  mbaDisplayConnectome(mbaFiberSplitLoops(fasAcpc.fibers),fig(3).h, [.1 .45 .95],'single',[],[],.2);
0150     view(-23,-23);delete(fig(3).light); fig(3).light = camlight('right');
0151         
0152     tmp_fas = feGet(feNoLesion,'fibers acpc');
0153     fibers2display = randsample(1:length(tmp_fas.fibers),ceil(0.01*length(tmp_fas.fibers)));
0154     fig(4).name = sprintf('pathneighborhood_%s',mfilename);
0155     fig(4).h   = figure('name',fig(4).name,'color',[.1 .45 .95]);     
0156     fig(4).type = 'jpg';
0157     [fig(4).h, fig(4).light] =  mbaDisplayConnectome(mbaFiberSplitLoops(...
0158         tmp_fas.fibers(fibers2display)),fig(4).h, [.95 .45 .1],'single',[],.6,.2);
0159     view(-23,-23);delete(fig(2).light); fig(4).light = camlight('right');
0160 
0161     % (2) The facicle with the pah neighborhood
0162     fig(5).name = sprintf('pathneighborhood_AND_tract_%s',mfilename);
0163     fig(5).h   = figure('name',fig(5).name,'color',[.1 .45 .95]);
0164     fig(5).type = 'jpg';
0165     [fig(5).h, fig(5).light] =  mbaDisplayConnectome(mbaFiberSplitLoops(fasAcpc.fibers),fig(5).h, [.1 .45 .95],'single',[],[],.2);
0166     view(-23,-23);delete(fig(5).light);
0167     hold on
0168     [fig(5).h, fig(5).light] =  mbaDisplayConnectome(mbaFiberSplitLoops(...
0169         tmp_fas.fibers(fibers2display)),fig(5).h, [.95 .45 .1],'single',[],.6,.2);
0170     view(-23,-23);delete(fig(5).light); fig(5).light = camlight('right');
0171  
0172     % Show the voxels coordinates to see if we are in the right spot
0173     plot_test_fiber_location_debug= false;
0174     if plot_test_fiber_location_debug
0175         figure('name','Coordinate check');
0176         plot3(allCoords(commonCoords,1),allCoords(commonCoords,2), ...
0177             allCoords(commonCoords,3),'ko','MarkerFaceColor','k','MarkerSize',8);
0178         hold on;
0179         plot3(fasCoords(:,1),fasCoords(:,2),fasCoords(:,3),'ro', ...
0180             'MarkerFaceColor',[.95 .45 .1],'MarkerSize',3);
0181         axis equal
0182         view(-23,-23); hold on
0183         plot3(allCoords(commonCoords,1),allCoords(commonCoords,2), ...
0184             allCoords(commonCoords,3),'go','MarkerFaceColor','g','MarkerSize',12);
0185         view(-23,-23);
0186     end
0187 end
0188 
0189 % RMSE distributions
0190 fig(6).name = sprintf('Size_of_effect_of_the_lesion_%s',mfilename);
0191 fig(6).type = 'eps';   
0192 fig(6).h   = nan;
0193 if display.evidence
0194     fig(6).name = sprintf('Size_of_effect_of_the_lesion_%s',mfilename);
0195     fig(6).h   = figure('name',fig(6).name,'color','w');
0196     set(fig(6).h,'Units','normalized','Position',[0.007 0.55  0.28 0.36]);
0197     subplot(1,4,1)
0198     plot(1,se.s.mean,'-o','color', [.95 .45 .1],'linewidth',2); hold on
0199     plot([1,1], [se.s.mean,se.s.mean] + [-se.s.std,se.s.std], '-','color',[.95 .45 .1] )
0200     ylabel('S (s.d.)', 'fontsize',14);
0201     set(gca,'box','off','xlim',[0 2], 'ylim',[0 ceil(se.s.mean + se.s.std)], ...
0202         'tickdir', 'out', 'ticklength', [0.025 0])
0203     subplot(1,4,2)
0204     plot(1,se.em.mean,'-o','color', [.95 .45 .1],'linewidth',2); hold on
0205     ylabel('Earth mover''s distance (raw scanner units)', 'fontsize',14);
0206     set(gca,'box','off','xlim',[0 2], 'ylim',[0 ceil(se.em.mean)], ...
0207         'tickdir', 'out', 'ticklength', [0.025 0])
0208     subplot(1,4,3)
0209     plot(1,se.kl.mean,'-o','color', [.95 .45 .1],'linewidth',2); hold on
0210     ylabel('K-L divergence (bits)', 'fontsize',14);
0211     set(gca,'box','off','xlim',[0 2], 'ylim',[0 ceil(se.kl.mean)], ...
0212         'tickdir', 'out', 'ticklength', [0.025 0])
0213     subplot(1,4,4)
0214     plot(1,se.j.mean,'-o','color', [.95 .45 .1],'linewidth',2); hold on
0215     ylabel('Jeffrey''s divergence (bits)', 'fontsize',14);
0216     set(gca,'box','off','xlim',[0 2], 'ylim',[0 ceil(se.j.mean)], ...
0217         'tickdir', 'out', 'ticklength', [0.025 0])
0218 end
0219 
0220 end