function mri_display(anat, tmap, scale, transp);

% function mri_display(anat, tmap, scale, transp);
%
% Aug 2003, Farshad Moradi
% California Institute of Technology
%
% To display anatomical images, use
%   mri_display(anat, 0, [0 3 8]);
%
% To superimpose a t-map for voxels with |t|>4,
%   mri_display(anat, tmap, [0 4 8]);
%
% If original scans are transverse, to view slices as sagital:
%   mri_display(anat, tmap, [0 4 8], 3);
% To view coronal slices:
%   mri_display(anat, tmap, [0 4 8], 4);
%

if nargin<4, transp = 1; end;
if nargin<3, scale = [0,3,8]; end;
if nargin<2, tmap = anat*0; end;

anat = conv3dToMosaic(anat, transp);
tmap = conv3dToMosaic(tmap, transp);

if size(anat,1)==2*size(tmap,1),
    tmap = tmap(floor(1:0.5:size(tmap,1)+0.5),:);
end;
if size(anat,2)==2*size(tmap,2),
    tmap = tmap(:,floor(1:0.5:size(tmap,2)+0.5));
end;

if tmap==0, tmap = anat*0; end;

if scale(1)<=0,
    if scale(1)<0, 
        anat = anat-min(anat(:)); 
        anat = 128*anat./max(1+anat(:));
    else,
        f = find(anat<0);
        tmap(f)=0;
        anat = 128*anat./max(1+anat(:));
        anat(f) = 127;
    end;
else,
    anat = 128*anat./scale(1);
    f = find(anat<0); tmap(f)=0; anat(f) = 127;
end;
f = find(abs(tmap)<=scale(2));
tmap(f)=0;
s = sign(tmap);
tmap = 64*(tmap - s*scale(2))/(scale(3)-scale(2));

ima = anat+65;
f = find(tmap>0); 
ima(f) = tmap(f)+193;
f = find(tmap<0);
ima(f) = 64+tmap(f);
cmap_t = jet(128+16);
cmap_a = gray(128);
clf;
if nargin>=2,
    SUBPLOT('position',[0.11 0.11 0.8 0.815]);
end;
f = find(ima<1); ima(f)=1;
image(ima); colormap([cmap_t(64:-1:1,:); cmap_a; cmap_t(end:-1:65+16,:)]); axis off;
if nargin>=2,
	SUBPLOT('position',[0.95 0.5 0.02 0.4]);
	image([1],[-4:8/(128+16):4],[[256:-1:193] ones(1,16)*128 [64:-1:0]]'); 
	set(gca, 'XTickLabel',[]);
	set(gca, 'XTick',[]);
	set(gca, 'YTick', [-4 -0.5 0.5 4]);
	set(gca, 'YTickLabel', [scale(3) scale(2) -scale(2) -scale(3)]);
    title('t');
end;

function temp = conv3dToMosaic(img, transp);

if size(img, 3)>1,
    if transp>1,
        img = shiftdim(img, floor(transp/2));
    end;
    [a,b, mn] = size(img);
    m = floor(sqrt(mn));
    n = ceil(mn / m);
    temp = zeros(a*m, b*n);
    for i=0:m-1,
        for j = 0:n-1;
            if i+j*m<mn,
                temp(i*a+(1:a), j*b+(1:b)) = img(:,:, i+j*m+1);
            end;
        end;
    end;
else,
    temp = img;
end;
if mod(transp,2) == 1,
    temp = temp';
end;
temp = temp(end:-1:1,:);