Description of fePlot

0001 function [uData, g] = fePlot(fe,plotType,varargin)
0002 % Gateway routine for plotting from the fascicle evaluation (fe) struct
0003 %
0004 %   [uData, g] = fePlot(fe,plotType,varargin)
0005 %
0006 % Example
0007 %   fePlot(fe,'quality of fit');
0008 %   fePlot(fe,'A matrix');
0009 %
0010 %   fePlot(fe,'bvecs');
0011 %
0012 %   fePlot(fe,'dsig measured');
0013 %
0014 % Copyright (2013-2014), Franco Pestilli, Stanford University, pestillifranco@gmail.com.
0015 %
0016 %-----------
0017 % Plot the fiber group in a matlab 3D mesh.
0018 % uData = fePlot(fe,'connectome')
0019 %-----------
0020 % Plots the quality of fit of the LiFE model
0021 % fData = ePlot(fe,'quality of fit');
0022 %-----------
0023 % Show an image of the LiFE model/matrix
0024 % uData = fePlot(fe,'model')
0025 %-----------
0026 % Plot the series of the measured diffusion signal reordered in such a way
0027 % that close angle in on the sphere are close together in the series. This
0028 % call uses the functionality in feGet(fe,dsig measured', varargin)
0029 % uData = fePlot(fe,'dsig measured')
0030 % uData = fePlot(fe,'dsig measured',voxelIndex)
0031 % uData = fePlot(fe,'dsig measured',coords)
0032 %-----------
0033 % Shows the bvecs on a sphere.
0034 % uData = dwiPlot(fe,'bvecs')
0035 %-----------
0036 % Plot a histogram of the length of the fibers in the connectome.
0037 % uData = fePlot(fe,'fiber length')
0038 %-----------
0039 %
0040 % Copyright (2013-2014), Franco Pestilli, Stanford University, pestillifranco@gmail.com.
0041 
0042 if notDefined('fe'),       error('''fe'' structure required.'); end
0043 if notDefined('plotType'), error('plotType required'); end
0044 
0045 % Open up your window for the current figure
0046 g = mrvNewGraphWin(sprintf('LiFE - %s - %s',upper(plotType),feGet(fe,'name')));
0047 uData = [];
0048 
0049 % Format the input parameter
0050 plotType = mrvParamFormat(plotType);
0051 
0052 % Find the plot
0053 switch plotType
0054   case {'2dhistogram'}
0055       % Makes a 2D density map.
0056       uData.x  = feGetRep(fe,varargin{1});
0057       uData.y  = feGetRep(fe,varargin{2});
0058       uData.ax.min   = min([uData.x,uData.y]);
0059       uData.ax.max   = max([uData.x,uData.y]);
0060       uData.ax.res   = 100;
0061       uData.ax.bins  = linspace(uData.ax.min,uData.ax.max,uData.ax.res);
0062       ymap = hist3([uData.x;uData.y]',{uData.ax.bins, uData.ax.bins});
0063       uData.sh = imagesc(flipud(ymap));
0064       colormap(flipud(hot)); 
0065       view(0,90);
0066       axis('square')
0067       
0068   case {'wmvolume'}
0069     % get the fiber density
0070     mmpVx = feGet(fe,'xform img 2 acpc');
0071     mmpVx = mmpVx(1:3,1:3);
0072     mmpVx = diag(mmpVx);
0073     vxNum = size(feGet(fe,'roi coords'),1);
0074     connectomeVol = vxNum*prod(mmpVx);% in mm^3
0075         
0076     % Total WM volume in the original ROI
0077     testVol = '/azure/scr1/frk/rois/life_right_occipital_example_large_no_cc_smaller.mat';
0078     ROI = dtiReadRoi(testVol);
0079     vxNum    = size(unique(floor(mrAnatXformCoords(feGet(fe,'xform acpc 2 img'),ROI.coords)),'rows'),1);
0080     totWMvol = vxNum*prod(mmpVx);% in mm^3
0081     uData = connectomeVol/totWMvol*100;
0082     close(g)
0083 
0084     g(1) = mrvNewGraphWin(sprintf('%s_percentWMVolume',feGet(fe,'name')));
0085     set(gcf,'color','w')
0086     
0087     % Volume of the ROI
0088     h = bar(uData,'r');
0089     set(gca,'ylim',[60 120],'xlim',[.5 1.5])
0090     ylabel('Percent white-matter volume')
0091     title(sprintf('Percent volume: %2.2f',uData))
0092     set(gca,  'ytick',[60 80 100 120], ...
0093       'box','off','tickDir','out')
0094     
0095     g(2) = mrvNewGraphWin(sprintf('%s_WMVolume',feGet(fe,'name')));
0096     set(gcf,'color','w')
0097     
0098     % Volume of the ROI
0099     h = bar([totWMvol, connectomeVol],'r');
0100     set(gca,'ylim',[0 46000],'xlim',[0 3])
0101     ylabel('White-matter volume (mm^3)')
0102     title(sprintf('Volume: total %6.0fmm^3, connectome %6.0fmm^3',totWMvol,connectomeVol))
0103     set(gca,  'ytick',[0 46000/2 46000], ...
0104         'xticklabel',{'Total WM','Connectome'}, ...
0105         'box','off','tickDir','out')
0106     
0107   case {'fiberdensitymap'}
0108     % Select a slice of brain in sagittal vie.
0109     % only Sagittal view is implemented
0110     if isempty(varargin),  slice = 30;
0111     else slice = varargin{1}; end
0112     if ~isempty(varargin) && length(varargin)==2
0113         fd = varargin{2};
0114     else
0115         % get the fiber density
0116         fd = (feGet(fe,'fiber density'));
0117     end
0118     
0119     close(g)
0120     g(1) = mrvNewGraphWin(sprintf('%s_FiberDensityMapFull',feGet(fe,'name')));
0121     set(gcf,'color','w')   
0122     img = feReplaceImageValues(nan(feGet(fe,'map size')),fd(:,1)',feGet(fe,'roiCoords'));
0123     maxfd = nanmax(img(:)); % This will be used tonormalize the fiber density plots
0124     surf(((fliplr(img(:,:,slice)./maxfd)'*255)),...
0125       'facecolor','texture','faceAlpha',1,'edgealpha',0,'CDataMapping','Direct');
0126     axis off; axis equal;
0127     set(gca,'ylim',[75 95],'xlim',[40 65])
0128     view(0,-90)
0129     
0130     cmap = colormap(jet(255));
0131     colorbar('ytick',[.25*size(cmap,1) .5*size(cmap,1) .75*size(cmap,1) size(cmap,1)],'yticklabel', ...
0132         {num2str(ceil(maxfd/8)) num2str(ceil(maxfd/4)) ...
0133         num2str(ceil(maxfd/2)) num2str(ceil(maxfd))},'tickdir','out')
0134     
0135     % Fiber density after life
0136     g(2) = mrvNewGraphWin(sprintf('%s_FiberDensityMapLiFE',feGet(fe,'name')));
0137     set(gcf,'color','w')
0138     img = feReplaceImageValues(nan(feGet(fe,'map size')),(fd(:,2))',feGet(fe,'roiCoords'));
0139     surf(fliplr(img(:,:,slice)./maxfd)'*255,...
0140         'facecolor','texture','faceAlpha',1,'edgealpha',0,'CDataMapping','Direct');
0141     axis off; axis equal;
0142     set(gca,'ylim',[75 95],'xlim',[40 65])
0143     % set(gca,'ylim',[100 130],'xlim',[50 85])
0144     
0145     %set(gca,'ylim',[70 100],'xlim',[40 65])
0146     view(0,-90)
0147     
0148     cmap = colormap(jet(255));
0149     colorbar('ytick',[.25*size(cmap,1) .5*size(cmap,1) .75*size(cmap,1) size(cmap,1)],'yticklabel', ...
0150         {num2str(ceil(maxfd/8)) num2str(ceil(maxfd/4)) ...
0151         num2str(ceil(maxfd/2)) num2str(ceil(maxfd))},'tickdir','out')
0152     
0153     % Weigth density (sum of weights)
0154     g(3) = mrvNewGraphWin(sprintf('%s_SumOfWeightsMapLiFE',feGet(fe,'name')));
0155     set(gcf,'color','w')
0156     img = feReplaceImageValues(nan(feGet(fe,'map size')),(fd(:,3))',feGet(fe,'roiCoords'));
0157     maxw = nanmax(img(:)); % This will be used tonormalize the fiber density plots
0158     minw = nanmin(img(:)); % This will be used tonormalize the fiber density plots
0159     surf(fliplr(img(:,:,slice+10)./maxw)'*255,...
0160         'facecolor','texture','faceAlpha',1,'edgealpha',0,'CDataMapping','Direct');
0161     axis off; axis equal;
0162     set(gca,'ylim',[75 95],'xlim',[40 65])
0163     %set(gca,'ylim',[70 100],'xlim',[40 65])
0164     %set(gca,'ylim',[100 130],'xlim',[50 85])
0165     view(0,-90)
0166     
0167     cmap = colormap(jet(255));
0168     colorbar('ytick',[0 .5*size(cmap,1) size(cmap,1)],'yticklabel', ...
0169       {num2str(minw)  num2str(maxw/2)  num2str(maxw)},'tickdir','out')
0170     
0171   case {'fiberdensityhist'}
0172     % get the fiber density
0173     fd = feGet(fe,'fiber density');
0174     close(g)
0175 
0176     g(1) = mrvNewGraphWin(sprintf('%s_FiberDensityHistFull_&_LiFE',feGet(fe,'name')));
0177     set(gcf,'color','w')
0178     
0179     % Fiber density before life
0180     edges = logspace(.5,3.2,100);
0181     linspace()
0182     centers = sqrt(edges(1:end-1).*edges(2:end));
0183     uData.full = histc(fd(:,1),edges)/size(fd,1)*100;
0184     h = bar(uData.full,'r');
0185     set(h,'edgecolor','r','linewidth',.01)
0186     set(get(h,'Children'),'FaceAlpha',.5,'EdgeAlpha',.5)
0187     
0188     % Fiber density after life
0189     hold on
0190     uData.life = histc(fd(:,2),edges)/size(fd,1)*100;
0191     h = bar(uData.life,'b');
0192     set(h,'edgecolor','b','linewidth',.01)
0193     set(get(h,'Children'),'FaceAlpha',.35,'EdgeAlpha',.35)
0194     set(gca,'ylim',[0 3],'xlim',[1 100])
0195     ticks = get(gca,'xtick'); 
0196     ylabel('Percent white-matter volume')
0197     xlabel('Number of fibers per voxel')
0198     set(gca,  'ytick',[0 1 2 3], ...
0199       'xticklabel', ceil(centers(ticks-1)) ,...
0200       'box','off','tickDir','out','xscale','lin')
0201     
0202     % Weigth density (sum of weights)
0203     g(2) = mrvNewGraphWin(sprintf('%s_SumOfWeightsHistLiFE',feGet(fe,'name')));
0204     set(gcf,'color','w')
0205     edges = logspace(-7,0,100);
0206     centers = sqrt(edges(1:end-1).*edges(2:end));
0207     y = histc(fd(:,3),edges)/size(fd,1)*100;
0208     h = bar(y,'k');
0209     set(h,'edgecolor','k','linewidth',.01)
0210     ylabel('Percent white-matter volume')
0211     xlabel('Sum of fascicles'' contribution to the voxel signal')
0212     set(gca,'ylim',[0 ceil(max(y))],'xlim',[.5 100])
0213     ticks = get(gca,'xtick');
0214     set(gca, 'ytick',[0 ceil(max(y))./2 ceil(max(y))], ...
0215       'xticklabel', ceil(1000000*centers(ticks-1))/1000000 ,...
0216       'box','off','tickDir','out','xscale','lin')
0217     
0218     % Weight density (mean of weights)
0219     g(3) = mrvNewGraphWin(sprintf('%s_MedianWeightsHistLiFE',feGet(fe,'name')));
0220     set(gcf,'color','w')
0221     edges = logspace(-7,0,100);
0222     centers = sqrt(edges(1:end-1).*edges(2:end));
0223     y = histc(fd(:,4),edges)/size(fd,1)*100;
0224     h = bar(y,'k');
0225     set(h,'edgecolor','k','linewidth',.01)
0226     ylabel('Percent white-matter volume')
0227     xlabel('Mean fascicles'' contribution to the voxel signal')
0228       set(gca,'ylim',[0 ceil(max(y))],'xlim',[.5 100])
0229     ticks = get(gca,'xtick');
0230     set(gca, 'ytick',[0 ceil(max(y))./2 ceil(max(y))], ...
0231       'xticklabel', ceil(1000000*centers(ticks-1))/1000000 ,...
0232       'box','off','tickDir','out','xscale','lin')
0233 
0234     % Weight density (var of weights)
0235     g(4) = mrvNewGraphWin(sprintf('%s_VarianceWeightsHistLiFE',feGet(fe,'name')));
0236     set(gcf,'color','w')
0237     edges = logspace(-7,0,100);
0238     centers = sqrt(edges(1:end-1).*edges(2:end));
0239     y = histc(fd(:,5),edges)/size(fd,1)*100;
0240     h = bar(y,'k');
0241     set(h,'edgecolor','k','linewidth',.01)
0242     ylabel('Percent white-matter volume')
0243     xlabel('Variance of fascicles'' contribution to the voxel signal')
0244     set(gca,'ylim',[0 ceil(max(y))],'xlim',[.5 100])
0245     ticks = get(gca,'xtick');
0246     set(gca, 'ytick',[0 ceil(max(y))./2 ceil(max(y))], ...
0247       'xticklabel', ceil(1000000*centers(ticks-1))/1000000 ,...
0248       'box','off','tickDir','out','xscale','lin')
0249     
0250   case {'rmseratiohistogram'}
0251     % Weigth density after life
0252     set(gcf,'color','w')
0253     R = (feGet(fe,'rmseratio'));
0254  
0255     edges = logspace(-.3,.6,100);
0256     centers = sqrt(edges(1:end-1).*edges(2:end));
0257     y = histc(R,edges)/length(R)*100;
0258     h = bar(y,'k');
0259     set(h,'edgecolor','k','linewidth',.01)
0260     ylabel('Percent white-matter volume')
0261     xlabel('R_{rmse}')
0262     set(gca,'ylim',[0 6],'xlim',[.5 100])
0263     ticks = get(gca,'xtick');
0264     set(gca, 'ytick',[0 3 6], ...
0265       'xticklabel', ceil(100*centers(ticks-1))/100 ,...
0266       'box','off','tickDir','out','xscale','lin')
0267     
0268   case {'map','maprep'}
0269     % Generate the requested map and save it to volume.
0270     if (length(plotType) > 3)
0271       % Use the xvalidation 'get' routine
0272       map = feValues2volume(feGetRep(fe, varargin{1}), ...
0273         feGet(fe,'roi coords'), ...
0274         feGetRep(fe,'map size'));
0275     else
0276       % use the default 'get' routine.
0277       map = feValues2volume(feGet(fe, varargin{1}), ...
0278         feGet(fe,'roi coords'), ...
0279         feGet(fe,'map size'));
0280     end
0281     
0282     % Get the slices to display (acpc) -120 to 120
0283     if (length(varargin)==1), slices = [-11:13];
0284     else slices = varargin{2};end
0285     
0286     % The xform is taken from the LiFE structure.
0287     mapXform = feGet(fe,'xform img2acpc');
0288     
0289     % Get the T1 anatomy from the session.
0290     ni         = niftiRead(feGet(fe,'anatomy file'));
0291     anatomyImg = double(ni.data);
0292     
0293     % Get the anatomy xform, from the file?
0294     anatomyXform = ni.qto_xyz;
0295     clear ni;
0296     
0297     % Get the range of the map to display from the data.
0298     clippingRange = [0 max(map(:))];%[max([0,min(map(:))]) max([max(map(:)),1])];
0299     mrAnatOverlayMontage(map, mapXform, anatomyImg, anatomyXform, [], clippingRange,slices);
0300     colormap('autumn');colorbar('location','eastoutside')
0301     title(upper(varargin{1}))
0302     
0303   case 'rmse'
0304     % Plot a histogram of the RMSE of the LiFE fit.
0305     % uData = fePlot(fe,'rmse')
0306     uData.totVarExp = feGet(fe,'totpve');
0307     uData.totrmse   = feGet(fe,'totalrmse');
0308     uData.rmsebyVox   = feGet(fe,'voxrmse');
0309     uData.rmsemedian = nanmedian(uData.rmsebyVox);
0310 
0311     % male a histogram of the length of the fibers
0312     [y,x] = hist(uData.rmsebyVox,ceil(length(uData.rmsebyVox)/10));
0313     bar(x,y,.65,'r')
0314     hold on
0315     plot([uData.rmsemedian uData.rmsemedian],[0 400],'k')
0316     title(sprintf('Total percent variance explained: %2.2f\nTotal RMSE: %2.2f',uData.totVarExp, uData.totrmse))
0317     xlabel('Root mean square error')
0318     ylabel('Number of occurrences')
0319      
0320   case 'pve'
0321     % Plot a histogram of the percent variance explained in each voxel with
0322     % the LiFE fit.
0323     % uData = fePlot(fe,'pve')
0324     uData.totVarExp = feGet(fe,'totpve');
0325     uData.totrmse   = feGet(fe,'totalrmse');
0326     uData.pve   = 100*feGet(fe,'voxr2zero');
0327     uData.pve( uData.pve < -1000 ) = -200; % if some of the values are -Inf set them to -200
0328     uData.medianpve = nanmedian(uData.pve);
0329 
0330     % male a histogram of the length of the fibers
0331     [y,x] = hist(uData.pve,ceil(length(uData.pve)/10));
0332     bar(x,y,1,'r')    
0333     hold on
0334     plot([uData.medianpve uData.medianpve],[0 400],'k')
0335     title(sprintf('Percent variance explained: %2.2f\nRMSE: %2.2f',uData.totVarExp, uData.totrmse))
0336     xlabel('Percent variance explained')
0337     ylabel('Number of occurrences')
0338 
0339   case 'rmserep'
0340     % Plot a histogram of the RMSE of the LiFE fit.
0341     % uData = fePlot(fe,'rmse')
0342     uData.totVarExp  = feGetRep(fe,'totpve');
0343     uData.totrmse    = feGetRep(fe,'totalrmse');
0344     uData.rmsebyVox  = feGetRep(fe,'voxrmse');
0345     uData.rmsemedian = nanmedian(uData.rmsebyVox);
0346     
0347     % male a histogram of the length of the fibers
0348     [y,x] = hist(uData.rmsebyVox,ceil(length(uData.rmsebyVox)/10));
0349     bar(x,y,.65,'r')
0350     hold on
0351     plot([uData.rmsemedian uData.rmsemedian],[0 400],'k')
0352     title(sprintf('Total percent variance explained: %2.2f\nTotal RMSE: %2.2f',uData.totVarExp, uData.totrmse))
0353     xlabel('Root mean square error')
0354     ylabel('Number of occurrences')
0355      
0356   case 'pverep'
0357     % Plot a histogram of the percent variance explained in each voxel with
0358     % the LiFE fit.
0359     % uData = fePlot(fe,'pve')
0360     uData.totVarExp = feGetRep(fe,'totpve');
0361     uData.totrmse   = feGetRep(fe,'totalrmse');
0362     uData.pve       = 100*feGetRep(fe,'voxr2zero');
0363     uData.pve( uData.pve < -1000 ) = -200; % if some of the values are -Inf set them to -200
0364     uData.medianpve = nanmedian(uData.pve);
0365 
0366     % male a histogram of the length of the fibers
0367     [y,x] = hist(uData.pve,ceil(length(uData.pve)/10));
0368     bar(x,y,1,'r')    
0369     hold on
0370     plot([uData.medianpve uData.medianpve],[0 400],'k')
0371     title(sprintf('Percent variance explained: %2.2f\nRMSE: %2.2f',uData.totVarExp, uData.totrmse))
0372     xlabel('Percent variance explained')
0373     ylabel('Number of occurrences')
0374      
0375   case 'rmserepdata'
0376     % Plot a histogram of the RMSE of the LiFE fit.
0377     % uData = fePlot(fe,'rmse')
0378     uData.totVarExp  = feGetRep(fe,'totpvedata');
0379     uData.totrmse    = feGetRep(fe,'totalrmsedata');
0380     uData.rmsebyVox  = feGetRep(fe,'voxrmsedata');
0381     uData.rmsemedian = nanmedian(uData.rmsebyVox);
0382      
0383     % male a histogram of the length of the fibers
0384     [y,x] = hist(uData.rmsebyVox,ceil(length(uData.rmsebyVox)/10));
0385     bar(x,y,.65,'g')
0386     hold on
0387     plot([uData.rmsemedian uData.rmsemedian],[0 400],'k')
0388     title(sprintf('Total percent variance explained by the data : %2.2f\nTotal RMSE: %2.2f',uData.totVarExp, uData.totrmse))
0389     xlabel('Root mean square error')
0390     ylabel('Number of occurrences')
0391   
0392   case 'pverepdata'
0393     % Plot a histogram of the percent variance explained in each voxel with
0394     % the LiFE fit.
0395     % uData = fePlot(fe,'pve')
0396     uData.totVarExp = feGetRep(fe,'totpvedata');
0397     uData.totrmse   = feGetRep(fe,'totalrmsedata');
0398     uData.pve       = 100*feGetRep(fe,'voxr2zerodata');
0399     uData.pve( uData.pve < -1000 ) = -200; % if some of the values are -Inf set them to -200
0400     uData.medianpve = nanmedian(uData.pve);
0401     
0402     % male a histogram of the length of the fibers
0403     [y,x] = hist(uData.pve,ceil(length(uData.pve)/10));
0404     bar(x,y,1,'g')    
0405     hold on
0406     plot([uData.medianpve uData.medianpve],[0 400],'k')
0407     title(sprintf('Percent variance explained by the data: %2.2f\nRMSE: %2.2f',uData.totVarExp, uData.totrmse))
0408     xlabel('Percent variance explained')
0409     ylabel('Number of occurrences')
0410   
0411   case 'rmseratio'
0412     % Plot a histogram of the RMSE of the LiFE fit.
0413     % uData = fePlot(fe,'rmseratio')
0414     uData.totVarExpD   = feGetRep(fe,'totpvedata');
0415     uData.totVarExpM  = feGetRep(fe,'totpve');
0416     uData.totVarExpMvw= 100*feGetRep(fe,'totalr2voxelwise');
0417     uData.totrmse     = feGetRep(fe,'totalrmseratio');
0418     uData.rmsebyVox   = feGetRep(fe,'voxrmseratio');
0419     uData.medianRmse  = nanmedian(uData.rmsebyVox);
0420     
0421     % male a histogram of the length of the fibers
0422     [y,x] = hist(uData.rmsebyVox,ceil(length(uData.rmsebyVox)/10));
0423     bar(x,y,.65,'k')
0424     hold on
0425     plot([uData.medianRmse uData.medianRmse],[0 400],'k')
0426     title(sprintf('PVE Data1/Data2 (%2.2f) | Global (%2.2f) | Voxel-wise (%2.2f)\nRMSE Ratio, Total (%2.2f) | Median (%2.2f)', ...
0427       uData.totVarExpD, ...
0428       uData.totVarExpM, ...
0429       uData.totVarExpMvw, ...
0430       uData.totrmse,   ...
0431       uData.medianRmse))
0432     xlabel('Root mean square error')
0433     ylabel('Number of occurrences')
0434 
0435   case 'rmseratiovoxelwise'
0436     % Plot a histogram of the RMSE of the LiFE fit.
0437     % uData = fePlot(fe,'rmseratio')
0438     uData.totVarExpD  = feGetRep(fe,'totpvedata');
0439     uData.totVarExpM  = feGetRep(fe,'totpve');
0440     uData.totVarExpMvw= 100*feGetRep(fe,'totalr2voxelwise');
0441     uData.totrmse     = feGetRep(fe,'totalrmseratiovoxelwise');
0442     uData.rmsebyVox   = feGetRep(fe,'voxrmseratiovoxelwise');
0443     uData.medianRmse  = nanmedian(uData.rmsebyVox);
0444     
0445     % male a histogram of the length of the fibers
0446     [y,x] = hist(uData.rmsebyVox,ceil(length(uData.rmsebyVox)/10));
0447     bar(x,y,.65,'k')
0448     hold on
0449     plot([uData.medianRmse uData.medianRmse],[0 400],'k')
0450     title(sprintf('PVE Data1/Data2 (%2.2f) | Global (%2.2f) | Voxel-wise (%2.2f)\nRMSE Ratio, Total (%2.2f) | Median (%2.2f)', ...
0451       uData.totVarExpD, ...
0452       uData.totVarExpM, ...
0453       uData.totVarExpMvw, ...
0454       uData.totrmse,   ...
0455       uData.medianRmse))
0456     xlabel('Root mean square error')
0457     ylabel('Number of occurrences')
0458 
0459   case 'fiberlength'
0460     % Plot a histogram of the length of the fibers in the connectome.
0461     % uData = fePlot(fe,'fiber length')
0462     uData.len       = fefgGet(feGet(fe,'fg acpc'),'length');
0463     uData.nanmedian = nanmedian(uData.len);
0464     
0465     % male a histogram of the length of the fibers
0466     [y,x] = hist(uData.len,ceil(length(uData.len)/10));
0467     bar(x,y,2,'r')
0468     hold on
0469     plot([uData.nanmedian uData.nanmedian],[0 1400],'k')
0470     title(sprintf('Total number of fibers: %i',feGet(fe,'n fibers')))
0471     xlabel('Fiber length in mm')
0472     ylabel('Number of occurrences')
0473     
0474   case 'qualityoffit'
0475     % Plots the quality of fit of the LiFE model
0476     % fData = fePlot(fe,'quality of fit');
0477     uData = fePlotQualityOfFit(fe);
0478     
0479   case {'mmatrix','model'}
0480     % Show an image of the LiFE model/matrix
0481     % uData = fePlot(fe,'model')
0482     uData.M = feGet(fe,'M fiber');
0483     imagesc(uData.M);
0484     
0485   case 'dsigmeasured'
0486     % Plot the series of the measured diffusion signal reordered in such a
0487     % way that close angle in on the sphere are close together in the
0488     % series. This call uses the functionality in feGet(fe,dsig measured',
0489     % varargin)
0490     % uData = fePlot(fe,'dsig measured')
0491     % uData = fePlot(fe,'dsig measured',voxelIndex)
0492     % uData = fePlot(fe,'dsig measured',coords)
0493     
0494     % Use routine in dwiPlot to reorder the x-axis of this so that
0495     % every voxel is plotted
0496     dsig = feGet(fe,'dsig measured');
0497     dsig = dsig(1:150);
0498     
0499     % From dwiPlot
0500     bvecs = feGet(fe,'bvecs');
0501     
0502     % Prepare the signal in 2D (imag) format:
0503     bFlat = sphere2flat(bvecs,'xy');
0504     
0505     % Interpolating function on the 2D representation
0506     F = TriScatteredInterp(bFlat,dsig(:));
0507     
0508     % Set the (x,y) range and interpolate
0509     x = linspace(min(bFlat(:,1)),max(bFlat(:,1)),sqrt(feGet(fe,'nbvecs')));
0510     y = linspace(min(bFlat(:,1)),max(bFlat(:,1)),sqrt(feGet(fe,'nbvecs')));
0511     [X Y] = meshgrid(x,y);
0512     est = F(X,Y);
0513     est(isnan(est)) = 0;  % Make extrapolated 0 rather than NaN
0514     %est = round(est./2);
0515     
0516     inData.x = X;
0517     inData.y = Y;
0518     inData.data = est;
0519     
0520     uData = dwiSpiralPlot(inData);
0521     dwiSpiralPlot(inData);
0522 
0523   case 'bvecs'
0524     % Shows the bvecs on a sphere.
0525     % uData = dwiPlot(fe,'bvecs')
0526     uData = dwiPlot(feGet(fe,'dwi'),'bvecs');
0527     
0528   case {'connectome','fg','fibers'}    
0529     % Plot the fiber group in a matlab 3D mesh.
0530     % uData = fePlot(fe,'fg')
0531     % We also look for fibers that are 'broken.' And display those as
0532     % separated fibers.
0533     uData = feConnectomeDisplay(feSplitLoopFibers( feGet(fe,'fg img')),g);
0534     %uData = feConnectomeDisplay(feGet(fe,'fg img'),figure);
0535 
0536   otherwise
0537     error('Unknown plot type %s\n',plotType);
0538 end
0539 
0540 % Set Figure properties
0541 if ~notDefined('uData')
0542   set(g,'userdata',uData)
0543 end
0544 
0545 % Set axes properties
0546 set(gca, 'FontSize', 16, ...
0547   'tickDir','out', ...
0548   'box','off');
0549 
0550 end
0551 
0552 
0553 %% Scratch
0554 % g = mrvNewGraphWin;
0555 % u = randn(128,128);
0556 % clear uData
0557 % uData.u = u;
0558 % imagesc(u);
0559 % set(g,'userdata',uData);
0560 %
0561 % foo = get(g,'userdata')
0562
fePlot

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
