0001 function [uData, g] = fePlot(fe,plotType,varargin)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042 if notDefined('fe'), error('''fe'' structure required.'); end
0043 if notDefined('plotType'), error('plotType required'); end
0044
0045
0046 g = mrvNewGraphWin(sprintf('LiFE - %s - %s',upper(plotType),feGet(fe,'name')));
0047 uData = [];
0048
0049
0050 plotType = mrvParamFormat(plotType);
0051
0052
0053 switch plotType
0054 case {'2dhistogram'}
0055
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
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);
0075
0076
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);
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
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
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
0109
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
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(:));
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
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
0144
0145
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
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(:));
0158 minw = nanmin(img(:));
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
0164
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
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
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
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
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
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
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
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
0270 if (length(plotType) > 3)
0271
0272 map = feValues2volume(feGetRep(fe, varargin{1}), ...
0273 feGet(fe,'roi coords'), ...
0274 feGetRep(fe,'map size'));
0275 else
0276
0277 map = feValues2volume(feGet(fe, varargin{1}), ...
0278 feGet(fe,'roi coords'), ...
0279 feGet(fe,'map size'));
0280 end
0281
0282
0283 if (length(varargin)==1), slices = [-11:13];
0284 else slices = varargin{2};end
0285
0286
0287 mapXform = feGet(fe,'xform img2acpc');
0288
0289
0290 ni = niftiRead(feGet(fe,'anatomy file'));
0291 anatomyImg = double(ni.data);
0292
0293
0294 anatomyXform = ni.qto_xyz;
0295 clear ni;
0296
0297
0298 clippingRange = [0 max(map(:))];
0299 mrAnatOverlayMontage(map, mapXform, anatomyImg, anatomyXform, [], clippingRange,slices);
0300 colormap('autumn');colorbar('location','eastoutside')
0301 title(upper(varargin{1}))
0302
0303 case 'rmse'
0304
0305
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
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
0322
0323
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;
0328 uData.medianpve = nanmedian(uData.pve);
0329
0330
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
0341
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
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
0358
0359
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;
0364 uData.medianpve = nanmedian(uData.pve);
0365
0366
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
0377
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
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
0394
0395
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;
0400 uData.medianpve = nanmedian(uData.pve);
0401
0402
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
0413
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
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
0437
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
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
0461
0462 uData.len = fefgGet(feGet(fe,'fg acpc'),'length');
0463 uData.nanmedian = nanmedian(uData.len);
0464
0465
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
0476
0477 uData = fePlotQualityOfFit(fe);
0478
0479 case {'mmatrix','model'}
0480
0481
0482 uData.M = feGet(fe,'M fiber');
0483 imagesc(uData.M);
0484
0485 case 'dsigmeasured'
0486
0487
0488
0489
0490
0491
0492
0493
0494
0495
0496 dsig = feGet(fe,'dsig measured');
0497 dsig = dsig(1:150);
0498
0499
0500 bvecs = feGet(fe,'bvecs');
0501
0502
0503 bFlat = sphere2flat(bvecs,'xy');
0504
0505
0506 F = TriScatteredInterp(bFlat,dsig(:));
0507
0508
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;
0514
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
0525
0526 uData = dwiPlot(feGet(fe,'dwi'),'bvecs');
0527
0528 case {'connectome','fg','fibers'}
0529
0530
0531
0532
0533 uData = feConnectomeDisplay(feSplitLoopFibers( feGet(fe,'fg img')),g);
0534
0535
0536 otherwise
0537 error('Unknown plot type %s\n',plotType);
0538 end
0539
0540
0541 if ~notDefined('uData')
0542 set(g,'userdata',uData)
0543 end
0544
0545
0546 set(gca, 'FontSize', 16, ...
0547 'tickDir','out', ...
0548 'box','off');
0549
0550 end
0551
0552
0553
0554
0555
0556
0557
0558
0559
0560
0561
0562