0001 function varargout = slice_overlay(action, 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
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082 global SO
0083
0084 if nargin < 1
0085 action = 'display';
0086 else
0087 action = lower(action);
0088 end
0089
0090 switch action
0091 case 'checkso'
0092 checkso;
0093 case 'getcmap'
0094 varargout = {getcmap(varargin{1})};
0095 case 'volmaxmin'
0096 [mx mn] = volmaxmin(varargin{1});
0097 varargout = {mx, mn};
0098 case 'display'
0099
0100
0101 checkso;
0102
0103
0104 X=1;Y=2;Z=3;
0105 dims = SO.slicedef;
0106 xmm = dims(X,1):dims(X,2):dims(X,3);
0107 ymm = dims(Y,1):dims(Y,2):dims(Y,3);
0108 zmm = SO.slices;
0109 [y x] = meshgrid(ymm,xmm');
0110 vdims = [length(xmm),length(ymm),length(zmm)];
0111
0112
0113 nslices = vdims(Z);
0114 minnpanels = nslices;
0115 cbars = 0;
0116 if is_there(SO,'cbar')
0117 cbars = length(SO.cbar);
0118 minnpanels = minnpanels+cbars;
0119 end
0120
0121
0122
0123 figno = figure(SO.figure);
0124
0125
0126
0127
0128 if ~SO.refreshf
0129 axisd = flipud(findobj(SO.figure, 'Type','axes','Tag', 'slice overlay panel'));
0130 npanels = length(axisd);
0131 if npanels < vdims(Z)+cbars;
0132 SO.refreshf = 1;
0133 end
0134 end
0135 if SO.refreshf
0136
0137 clf
0138 axisd = [];
0139
0140
0141 set(figno,'InvertHardCopy','off');
0142
0143
0144 parea = SO.area.position;
0145 if ~strcmp(SO.area.units, 'pixels')
0146 ubu = get(SO.figure, 'units');
0147 set(SO.figure, 'units','pixels');
0148 tmp = get(SO.figure, 'Position');
0149 ascf = tmp(3:4);
0150 if ~strcmp(SO.area.units, 'normalized')
0151 set(SO.figure, 'units',SO.area.units);
0152 tmp = get(SO.figure, 'Position');
0153 ascf = ascf ./ tmp(3:4);
0154 end
0155 set(figno, 'Units', ubu);
0156 parea = parea .* repmat(ascf, 1, 2);
0157 end
0158 asz = parea(3:4);
0159
0160
0161 yxratio = length(ymm)*dims(Y,2)/(length(xmm)*dims(X,2));
0162 if ~is_there(SO, 'xslices')
0163
0164 axlen(X,:)=asz(1):-1:1;
0165 axlen(Y,:)=yxratio*axlen(X,:);
0166 panels = floor(asz'*ones(1,size(axlen,2))./axlen);
0167 estnpanels = prod(panels);
0168 tmp = find(estnpanels >= minnpanels);
0169 if isempty(tmp)
0170 error('Whoops, cannot fit panels onto figure');
0171 end
0172 b = tmp(1);
0173 panels = panels(:,b);
0174 axlen = axlen(:, b);
0175 else
0176
0177 panels(X) = SO.xsliceno;
0178 axlen(X) = asz(X)/panels(X);
0179 end
0180
0181
0182 panels(Y) = ceil(minnpanels/panels(X));
0183 axlen(Y) = axlen(X)*yxratio;
0184
0185
0186 divs = [Inf 2 1];the_ds = [0;0];
0187 the_ds(X) = divs(strcmp(SO.area.halign, {'left','center','right'}));
0188 the_ds(Y) = divs(strcmp(SO.area.valign, {'bottom','middle','top'}));
0189 startc = parea(1:2)' + (asz'-(axlen.*panels))./the_ds;
0190
0191
0192 r=0;c=1;
0193 npanels = prod(panels);
0194 lastempty = npanels-cbars;
0195 for i = 1:npanels
0196
0197 if i<=nslices
0198 u.type = 'slice';
0199 u.no = zmm(i);
0200 elseif i > lastempty
0201 u.type = 'cbar';
0202 u.no = i - lastempty;
0203 else
0204 u.type = 'empty';
0205 u.no = i - nslices;
0206 end
0207 axpos = [r*axlen(X)+startc(X) (panels(Y)-c)*axlen(Y)+startc(Y) axlen'];
0208 axisd(i) = axes(...
0209 'Parent',figno,...
0210 'XTick',[],...
0211 'XTickLabel',[],...
0212 'YTick',[],...
0213 'YTickLabel',[],...
0214 'Box','on',...
0215 'XLim',[1 vdims(X)],...
0216 'YLim',[1 vdims(Y)],...
0217 'Units', 'pixels',...
0218 'Position',axpos,...
0219 'Tag','slice overlay panel',...
0220 'UserData',u);
0221 r = r+1;
0222 if r >= panels(X)
0223 r = 0;
0224 c = c+1;
0225 end
0226 end
0227 end
0228
0229
0230 if is_there(SO,'labels')
0231 labels = SO.labels;
0232 if iscell(labels.format)
0233 if length(labels.format)~=vdims(Z)
0234 error(...
0235 sprintf('Oh dear, expecting %d labels, but found %d',...
0236 vdims(Z), length(labels.contents)));
0237 end
0238 else
0239
0240 fstr = labels.format;
0241 labels.format = cell(vdims(Z),1);
0242 acpt = SO.transform * [0 0 0 1]';
0243 for i = 1:vdims(Z)
0244 labels.format(i) = {sprintf(fstr,zmm(i)-acpt(Z))};
0245 end
0246 end
0247 end
0248
0249
0250 nimgs = length(SO.img);
0251 lrn = zeros(nimgs,3);
0252 cmaps = cell(nimgs);
0253 for i = 1:nimgs
0254 cmaps(i)={SO.img(i).cmap};
0255 lrnv = {SO.img(i).outofrange{:}, SO.img(i).nancol};
0256 for j = 1:length(lrnv)
0257 if prod(size(lrnv{j}))==1
0258 lrn(i,j) = lrnv{j};
0259 else
0260 cmaps(i) = {[cmaps{i}; lrnv{j}(1:3)]};
0261 lrn(i,j) = size(cmaps{i},1);
0262 end
0263 end
0264 end
0265
0266
0267 nvox = prod(vdims(1:2));
0268 pandims = [vdims([2 1]) 3];
0269
0270 zimg = zeros(pandims);
0271 for i = 1:nslices
0272 ixyzmm = [x(:)';y(:)';ones(1,nvox)*zmm(i);ones(1,nvox)];
0273 img = zimg;
0274
0275
0276 for j = 1:nimgs
0277 if ~any(ismember(j,SO.contours))
0278 thisimg = SO.img(j);
0279
0280
0281 vixyz = inv(SO.transform*thisimg.vol.mat)*ixyzmm;
0282
0283 i1 = spm_sample_vol(thisimg.vol,vixyz(X,:),vixyz(Y,:),vixyz(Z,:), ...
0284 [thisimg.hold thisimg.background]);
0285 if is_there(thisimg, 'func')
0286 eval(thisimg.func);
0287 end
0288
0289 i1 = reshape(i1, vdims(1:2))';
0290
0291 [csdata badvals]= scaletocmap(...
0292 i1,...
0293 thisimg.range(1),...
0294 thisimg.range(2),...
0295 cmaps{j},...
0296 lrn(j,:));
0297
0298 iimg = reshape(cmaps{j}(csdata(:),:),pandims);
0299 tmp = repmat(logical(~badvals),[1 1 3]);
0300 if thisimg.prop ~= Inf
0301 img(tmp) = img(tmp) + iimg(tmp)*thisimg.prop;
0302 else
0303 img(tmp) = iimg(tmp);
0304 end
0305 end
0306 end
0307
0308
0309 img(img>1) = 1;
0310
0311 image('Parent', axisd(i),...
0312 'CData',img);
0313 if is_there(SO,'labels')
0314 text('Parent',axisd(i),...
0315 'Color', labels.colour,...
0316 'FontUnits', 'normalized',...
0317 'VerticalAlignment','bottom',...
0318 'HorizontalAlignment','left',...
0319 'Position', [1 1],...
0320 'FontSize',labels.size,...
0321 'String', labels.format{i});
0322 end
0323
0324
0325 for j = 1:nimgs
0326 if any(ismember(j,SO.contours))
0327 thisimg = SO.img(j);
0328
0329 vixyz = inv(SO.transform*thisimg.vol.mat)*ixyzmm;
0330
0331 i1 = spm_sample_vol(thisimg.vol,vixyz(X,:),vixyz(Y,:),vixyz(Z,:), ...
0332 [thisimg.hold thisimg.background]);
0333 if is_there(thisimg, 'func')
0334 eval(thisimg.func);
0335 end
0336
0337 i1 = reshape(i1, vdims(1:2))';
0338
0339
0340 axes(axisd(i));
0341 set(axisd(i),'NextPlot','add');
0342 contour(i1,[thisimg.range(1) thisimg.range(1)],thisimg.linestr)
0343 end
0344 end
0345 end
0346
0347 for i = (nslices+1):npanels
0348 set(axisd(i),'Color',[0 0 0]);
0349 end
0350
0351 for i = 1:cbars
0352 axno = axisd(end-cbars+i);
0353 cbari = SO.img(SO.cbar(i));
0354 cml = size(cbari.cmap,1);
0355 p = get(axno, 'Position');;
0356 cw = p(3)*0.2;
0357 ch = p(4)*0.75;
0358 pc = p(3:4)/2;
0359 [axlims idxs] = sort(cbari.range);
0360 a=axes(...
0361 'Parent',figno,...
0362 'XTick',[],...
0363 'XTickLabel',[],...
0364 'Units', 'pixels',...
0365 'YLim', axlims,...
0366 'FontUnits', 'normalized',...
0367 'FontSize', 0.075,...
0368 'YColor',[1 1 1],...
0369 'Tag', 'cbar',...
0370 'Box', 'off',...
0371 'Position',[p(1)+pc(1)-cw/2,p(2)+pc(2)-ch/2,cw,ch]...
0372 );
0373 ih = image('Parent', a,...
0374 'YData', axlims(idxs),...
0375 'CData', reshape(cbari.cmap,[cml,1,3]));
0376
0377 end
0378
0379
0380 end
0381
0382 return
0383
0384 function checkso
0385
0386 global SO
0387
0388
0389 if is_there(SO, 'figure')
0390 try
0391 if ~strcmp(get(SO.figure,'Type'),'figure')
0392 error('Figure handle is not a figure')
0393 end
0394 catch
0395 error('Figure handle is not a valid figure')
0396 end
0397 else
0398
0399 SO.figure = spm_figure('FindWin', 'Graphics');
0400 if isempty(SO.figure)
0401 SO.figure = gcf;
0402 end
0403 end
0404
0405 if strcmp(get(SO.figure, 'Tag'),'Graphics')
0406 SO.area.position = [0 0 1 0.92];
0407 SO.area.units = 'normalized';
0408 SO.area.valign = 'top';
0409 end
0410
0411
0412 orientn = [];
0413 SO = set_def(SO, 'transform', 'axial');
0414 if ischar(SO.transform)
0415 orientn = find(strcmpi(SO.transform, {'axial','coronal','sagittal'}));
0416 if isempty(orientn)
0417 error(sprintf('Unexpected orientation %s', SO.transform));
0418 end
0419 ts = [0 0 0 0 0 0 1 1 1;...
0420 0 0 0 pi/2 0 0 1 -1 1;...
0421 0 0 0 pi/2 0 -pi/2 -1 1 1];
0422 SO.transform = spm_matrix(ts(orientn,:));
0423 end
0424
0425 if ~is_there(SO,'slicedef' | ~is_there(SO, 'slices'))
0426
0427 V = SO.img(1).vol;
0428 D = V.dim(1:3);
0429 T = SO.transform * V.mat;
0430 vcorners = [1 1 1; D(1) 1 1; 1 D(2) 1; D(1:2) 1; ...
0431 1 1 D(3); D(1) 1 D(3); 1 D(2:3) ; D(1:3)]';
0432 corners = T * [vcorners; ones(1,8)];
0433 SC = sort(corners');
0434 vxsz = sqrt(sum(T(1:3,1:3).^2));
0435
0436 SO = set_def(SO, 'slicedef',...
0437 [SC(1,1) vxsz(1) SC(8,1);SC(1,2) vxsz(2) SC(8,2)]);
0438 SO = set_def(SO, 'slices',[SC(1,3):vxsz(3):SC(8,3)]);
0439 end
0440
0441
0442 SO = set_def(SO, 'cbars', []);
0443
0444
0445 SO = set_def(SO, 'refreshf', 1);
0446
0447
0448 defstruct = struct('colour',[1 1 1],'size',0.075,'format', '%+30f');
0449 if ~isfield(SO, 'labels')
0450 SO.labels = defstruct;
0451 elseif ~isempty(SO.labels)
0452
0453 SO.labels = set_def(SO.labels, 'colour', defstruct.colour);
0454
0455 SO.labels = set_def(SO.labels, 'size', defstruct.size);
0456
0457 SO.labels = set_def(SO.labels, 'format', defstruct.format);
0458 end
0459
0460
0461 defarea = struct('position',[0 0 1 1],'units','normalized');
0462 SO = set_def(SO, 'area', defarea);
0463 if ~is_there(SO.area, 'position')
0464 SO.area = defarea;
0465 end
0466 if ~is_there(SO.area,'units')
0467 if (all(SO.area.position>=0 & SO.area.position<=1))
0468 SO.area.units = 'normalized';
0469 else
0470 SO.area.units = 'pixels';
0471 end
0472 end
0473 SO.area = set_def(SO.area,'halign', 'center');
0474 SO.area = set_def(SO.area,'valign', 'middle');
0475
0476
0477
0478
0479
0480 remcol = 1;
0481 for i = 1:length(SO.img)
0482 if ~is_there(SO.img(i),'hold')
0483 SO.img(i).hold = 1;
0484 end
0485 if ~is_there(SO.img(i),'background')
0486 SO.img(i).background = NaN;
0487 end
0488 if ~is_there(SO.img(i),'prop')
0489
0490 SO.img(i).prop = remcol/(length(SO.img)-i+1);
0491 remcol = remcol - SO.img(i).prop;
0492 end
0493 if ~is_there(SO.img(i),'range')
0494 [mx mn] = volmaxmin(SO.img(i).vol);
0495 SO.img(i).range = [mn mx];
0496 end
0497 if ~is_there(SO.img(i),'cmap')
0498 if SO.img(i).prop == Inf;
0499 SO.img(i).cmap = check_map(i, 'spm_cols');
0500 else
0501 SO.img(i).cmap = check_map(i, 'actc');
0502 end
0503 end
0504 if ~is_there(SO.img(i),'outofrange')
0505
0506 if SO.img(i).prop == Inf
0507 if xor(SO.img(i).range(1) < SO.img(i).range(2), ...
0508 SO.img(i).range(2) < 0)
0509 SO.img(i).outofrange = {[0],size(SO.img(i).cmap,1)};
0510 else
0511 SO.img(imgno).outofrange={[1], [0]};
0512 end
0513 else
0514 SO.img(i).outofrange = {1,size(SO.img(i).cmap,1)};
0515 end
0516 end
0517 for j=1:2
0518 if isempty(SO.img(i).outofrange{j})
0519 SO.img(i).outofrange(j) = {0};
0520 end
0521 end
0522 if ~is_there(SO.img(i),'nancol')
0523 SO.img(i).nancol = 0;
0524 end
0525 end
0526 return
0527
0528 function tf = is_there(a, fname)
0529
0530 tf = isfield(a, fname);
0531 if tf
0532 tf = ~isempty(getfield(a, fname));
0533 end
0534 return
0535
0536 function [img, badvals]=scaletocmap(inpimg,mn,mx,cmap,lrn)
0537 img = (inpimg-mn)/(mx-mn);
0538 cml = size(cmap,1);
0539 if cml==1
0540 img(img>=0 & img<=1)=1;
0541 else
0542 img = img*(cml-1)+1;
0543 end
0544 outvals = {img<1, img>cml, isnan(img)};
0545 img= round(img);
0546 badvals = zeros(size(img));
0547 for i = 1:length(lrn)
0548 if lrn(i)
0549 img(outvals{i}) = lrn(i);
0550 else
0551 badvals = badvals | outvals{i};
0552 img(outvals{i}) = 1;
0553 end
0554 end
0555 return
0556
0557 function st = set_def(st, fld, def)
0558 if ~is_there(st, fld)
0559 st = setfield(st, fld, def);
0560 end
0561
0562 function cm = check_map(idx, mapstr)
0563 try
0564 load(mapstr)
0565 catch
0566 error(sprintf(['No colormap defined for image %d, '...
0567 'and default map %s is not on the path'],...
0568 idx, mapstr));
0569 end
0570 cm = eval(mapstr);
0571 return
0572
0573 function [mx,mn] = volmaxmin(vol)
0574 mx = -Inf;mn=Inf;
0575 for i=1:vol.dim(3),
0576 tmp = spm_slice_vol(vol,spm_matrix([0 0 i]),vol.dim(1:2),[0 NaN]);
0577 tmp = tmp(find(finite(tmp(:))));
0578 if ~isempty(tmp)
0579 mx = max([mx; tmp]);
0580 mn = min([mn; tmp]);
0581 end
0582 end
0583 return
0584
0585 function cmap = getcmap(acmapname)
0586
0587 if ~isempty(acmapname)
0588 cmap = evalin('base',acmapname,'[]');
0589 if isempty(cmap)
0590
0591 tmp = strcmp(acmapname, {'red','green','blue'});
0592 if any(tmp)
0593 cmap = zeros(64,3);
0594 cmap(:,tmp) = ((0:63)/63)';
0595 else
0596
0597 [p f e] = fileparts(acmapname);
0598 acmat = fullfile(p, [f '.mat']);
0599 if exist(acmat, 'file')
0600 s = struct2cell(load(acmat));
0601 cmap = s{1};
0602 end
0603 end
0604 end
0605 end
0606 if size(cmap, 2)~=3
0607 warning('Colormap was not an N by 3 matrix')
0608 cmap = [];
0609 end