mfbox_plotic.m

toolbox for spm 5 for data, model free analysis

function fh=mfbox_plotic(component,timecourse,ictitle,stdfact,bg,atlas,outputopts)
% fh = mfbox_icplot(component,timecourse,ictitle,stdfact,bg,atlas,outputopts)
%
% component   = 3d matrix der komponente
%               (is cell contains component, spatial error, and temporal error)
% bg          = 3d matrix des hirnhintergrundes
% stdfact     = factor of the standard deviation to consider voxels as background noise
%               (if cell contains stdfact to hide from plot and remmean)
% atlas       = plot atlas (if cell contains fron,side and top image)
%
% fh = if output generate SVG of plot
%
% Copyright by Peter Gruber, Fabian J. Theis and Ingo R. Keck
% Signal Processing & Information Theory group
% Institute of Biophysics, University of Regensburg, Germany
% Homepage: http://research.fabian.theis.name
%           http://www-aglang.uni-regensburg.de
%
% This file is free software, subject to the
% GNU GENERAL PUBLIC LICENSE, see gpl.txt

if (nargin<7), outputopts = 0; end
if (nargin<6), atlas = 0; end

outputsvg = 0;
rmean = true;
design = {};

if (iscell(atlas))
    front_atlas = atlas{1};
    side_atlas = atlas{2};
    top_atlas = atlas{3};
    atlas = 0.5;
end
if (iscell(timecourse))
    if (length(timecourse)>2), design = timecourse(3:end); end
    timeline = timecourse{2};
    timecourse = timecourse{1};
else
    timeline = 1:length(timecourse);
end
if (any(size(timecourse)==1)), timecourse = timecourse(:); end

fighandle = 0;
if (nargout>0)
    if (iscell(outputopts))
        fighandle = outputopts{1};
        outputopts = outputopts{2};
    elseif (outputopts<2)
        outputsvg = outputopts+1;
    end
end

red = [1;0;0];
green = [0;1;0];
blue = [0;0;1];
white = [1;1;1];
black = [0;0;0];

if (iscell(stdfact)), rmean = stdfact{2}; stdfact = stdfact{1}; end

if (iscell(component) && length(component)==3)
    snv = component{2}; tnv = component{3};
    component = component{1};
    v = ~isnan(component);
    snv(snv~=0) = 1-(1+(abs(component(snv~=0))./snv(snv~=0))).^(-1);
    snv(snv==0) = 1;
elseif (iscell(component) && length(component)==2)
    snv = component{2};
    tnv = [];
    component = component{1};
    v = ~isnan(component);
    snv(snv~=0) = 1-(1+(abs(component(snv~=0))./snv(snv~=0))).^(-1);
    snv(snv==0) = 1;
elseif (iscell(component))
    component = component{1};
    v = ~isnan(component);
    tnv = [];
    snv = abs(var(component(v)));
    snv = 1-(1+(abs(component)/snv)).^(-1);
else
    v = ~isnan(component);
    tnv = [];
    snv = abs(var(component(v)));
    snv = 1-(1+(abs(component)/snv)).^(-1);
end

if (rmean), component(v) = component(v)-mean(component(v)); end
minic = min(component(v));
maxic = max(component(v));
minic = -max(abs(minic),abs(maxic));
maxic = -minic;

if (stdfact>0)
    va = var(component(v));
    component(v & ((component.^2)<va*(stdfact.^2))) = 0;
end
component(~v) = 0;
component = component/maxic; %/(1-atlas);

if (~isempty(bg))
    minbg = min(bg(:));
    maxbg = max(bg(:));
    if (maxbg>minbg)
        bg = (bg-minbg)/(maxbg-minbg)/3;
    end
end

if (outputsvg>0)
    fh = {};
else
    if (fighandle~=0)
        fh = fighandle;
        figure(fh);
        clf;
        set(fh,'Visible','Off');
    else
        fh = figure('Color',white,'Name',ictitle,'Visible','off');
    end
end

pdim = 2; %front
[c,t,ec] = mfbox_argmax(abs(component),pdim,component,snv);
c = c.*sign(t);
if (atlas>0)
    imgdim = size(front_atlas'); %182,182
else
    imgdim = size(c);
end
c = c';
ec = ec';
if (~isempty(bg))
    b = scaleimg(squeeze(max(bg,[],pdim)),imgdim,false)';
    if (atlas>0)
        bb = min(1-atlas*front_atlas,1-(1-atlas)*b);
    else
        bb = 1-b;
    end
else
    bb = [];
end
if (outputsvg>0)
    fh{end+1} = sprintf('<g transform="translate(43,46)">\n<g transform="rotate(180)">\n%s\n</g>\n<g transform="translate(-20,-41)">\n<text x="0" y="0" font-family="Verdana" font-size="4" text-anchor="middle">front</text>\n</g>\n</g>\n', ...
        makesvgimg(40,40,bb,c,ec,red,blue));
else
    bb = makepixelimg(bb,c,ec,red,blue);
    front = axes('Position',[0.03,0.54,0.4,0.4]);
    image(min(max(bb,0),1),'Parent',front);
    set(front,'YDir','normal','XDir','reverse','XTickMode','manual', ...
        'XTickLabelMode','manual','YTickMode','manual', ...
        'YTickLabelMode','manual','XTick',[],'XTickLabel',{}, ...
        'YTick',[],'YTickLabel',{});
    set(get(front,'Title'),'String',{'front'});
end

pdim = 1; %side
[c,t,ec] = mfbox_argmax(abs(component),pdim,component,snv);
c = c.*sign(t);
if (atlas>0)
    imgdim = size(side_atlas'); %218,182
else
    imgdim = size(c);
end
c = c';
ec = ec';
if (~isempty(bg))
    b = scaleimg(squeeze(max(bg,[],pdim)),imgdim,false)';
    if (atlas>0)
        bb = min(1-atlas*side_atlas,1-(1-atlas)*b);
    else
        bb = 1-b;
    end
else
    bb = [];
end
if (outputsvg>0)
    fh{end+1} = sprintf('<g transform="translate(98,46)">\n<g transform="rotate(180)">\n%s\n</g>\n<g transform="translate(-24,-41)">\n<text x="0" y="0" font-family="Verdana" font-size="4" text-anchor="middle">side</text>\n</g>\n</g>\n', ...
        makesvgimg(48,40,bb,c,ec,red,blue));
else
    bb = makepixelimg(bb,c,ec,red,blue);
    side = axes('Position',[0.5,0.54,0.48,0.4]);
    image(min(max(bb,0),1),'Parent',side);
    set(side,'YDir','normal','XDir','reverse','XTickMode','manual', ...
        'XTickLabelMode','manual','YTickMode','manual', ...
        'YTickLabelMode','manual','XTick',[],'XTickLabel',{}, ...
        'YTick',[],'YTickLabel',{});
    set(get(side,'Title'),'String',{'side'});
end

pdim = 3; %top
[c,t,ec] = mfbox_argmax(abs(component),pdim,component,snv);
c = c.*sign(t);
if (atlas>0)
    imgdim = size(top_atlas'); %182,218
else
    imgdim = size(c);
end
c = c';
ec = ec';
if (~isempty(bg))
    b = scaleimg(squeeze(max(bg,[],pdim)),imgdim,false)';
    if (atlas>0)
        bb = min(1-atlas*top_atlas,1-(1-atlas)*b);
    else
        bb = 1-b;
    end
else
    bb = [];
end
if (outputsvg>0)
    fh{end+1} = sprintf('<g transform="translate(43,99)">\n<g transform="rotate(180)">\n%s"\n</g>\n<g transform="translate(-20,-49)">\n<text x="0" y="0" font-family="Verdana" font-size="4" text-anchor="middle">top</text>\n</g>\n</g>\n', ...
        makesvgimg(40,48,bb,c,ec,red,blue));
else
    bb = makepixelimg(bb,c,ec,red,blue);
    top = axes('Position',[0.03,0.01,0.4,0.48]);
    image(min(max(bb,0),1),'Parent',top);
    set(top,'YDir','normal','XDir','reverse','XTickMode','manual', ...
        'XTickLabelMode','manual','YTickMode','manual', ...
        'YTickLabelMode','manual','XTick',[],'XTickLabel',{}, ...
        'YTick',[],'YTickLabel',{});
    set(get(top,'Title'),'String',{'top'});
end

scale = reshape([(1-blue)*(0:0.05:1)+blue*ones(1,21),(1-red)*(1:-0.05:0)+ ...
    red*ones(1,21)]',[1,42,3]);
if (outputsvg>0)
    ctmp = tempname();
    imwrite(scale,ctmp,'png');
    fh{end+1} = sprintf('<g transform="translate(51,47.5)">\n<image x="10" y="0" width="29" height="2.5" xlink:href="data:image/png;base64,%s"/>\n<g transform="translate(9,2)">\n<text x="0" y="0" font-size="2" font-family="Verdana" text-anchor="end">%1.1f</text>\n</g>\n<g transform="translate(40,2)">\n<text x="0" y="0" font-size="2" font-family="Verdana" text-anchor="start">%1.1f</text>\n</g>\n</g>', ...
        mfbox_base64([],ctmp),minic,maxic);
    delete(ctmp);
else
    colorc = axes('Position',[0.58,0.425,0.35,0.025]);
    image(scale,'Parent',colorc);
    set(colorc,'XTickMode','manual', ...
        'XTickLabelMode','manual','YTickMode','manual', ...
        'YTickLabelMode','manual','YTick',[],'YTickLabel',{}, ...
        'XTick',[1;42],'XTickLabel',{sprintf('%1.1f',minic) sprintf('%1.1f',maxic)});
end

bb = []; c = ones(1,11); ec = [0:0.1:1];
if (outputsvg>0)
    fh{end+1} = sprintf('<g transform="translate(51,52.5)">\n<g transform="translate(10,0)">%s</g>\n<g transform="translate(9,2)">\n<text x="0" y="0" font-size="2" font-family="Verdana" text-anchor="end">%1.1f</text>\n</g>\n<g transform="translate(40,2)">\n<text x="0" y="0" font-size="2" font-family="Verdana" text-anchor="start">%1.1f</text>\n</g>\n</g>\n', ...
        makesvgimg(29,2.5,bb,c,ec,black,black),0,1);
else
    bb = makepixelimg(bb,c,ec,black,black);
    top = axes('Position',[0.58,0.5,0.35,0.025]);
    image(min(max(bb,0),1),'Parent',top);
    set(top,'XTickMode','manual', ...
        'XTickLabelMode','manual','YTickMode','manual', ...
        'YTickLabelMode','manual','YTick',[],'YTickLabel',{}, ...
        'XTick',[1;42],'XTickLabel',{sprintf('%1.1f',0) sprintf('%1.1f',1)});
end

mi = min(timecourse(:));
ma = max(timecourse(:));
if (~isempty(design))
    t = reshape(cat(3,design{:}),[],length(design))>0.5;
    c = repmat(reshape([red,green,blue]',[3,1,3]),[ceil(size(t,2)/3),size(t,1),1]);
    c = c(1:size(t,2),:,:);
    c = c.*repmat(reshape(t',[size(t,2),size(t,1),1]),[1,1,3]);
    cc = 1-0.25*c;
else
    cc = reshape(repmat([0.7;0.95],size(timecourse,1),1)*(red+green+blue)',1,2*size(timecourse,1),3);
end
if (outputsvg>0)
    atmp = tempname();
    imwrite(max(min(cc,1),0),atmp,'png');
    fh{end+1} = sprintf('<g transform="translate(50,70)">\n<rect x="5" y="0" width="42.5" height="25" fill="none" stroke="black" stroke-width="0.5"/>\n<image x="%1.4f" y="5" width="%1.4f" height="17" xlink:href="data:image/png;base64,%s"/>', ...
        7.5-1/size(cc,2),37.5*(size(cc,2)+2)/size(cc,2),mfbox_base64([],atmp));
    delete(atmp);
    if (ma>mi)
        tl = (timeline(:)-timeline(1))/(timeline(end)-timeline(1));
        for i=1:size(timecourse,2)
            tc = (timecourse(:,i)-mi)/(ma-mi);
            v = [37.5*tl,25*(1-tc)]';
            fh{end+1} = sprintf('<g transform="translate(7.5,1)">\n<polyline fill="none" stroke="red" stroke-width="0.5" points="%s"/>\n</g>', ...
                sprintf('%1.4f,%1.4f ',v(:)));
        end
        v = linspace(mi,ma,4); tc = linspace(26,1,4);
        for i=1:4
            fh{end+1} = sprintf('<line x1="5" x2="6" y1="%1.4f" y2="%1.4f" stroke="black" stroke-width="0.2"/>\n<g transform="translate(4,%1.4f)">\n<text x="0" y="0" font-size="2" font-family="Verdana" text-anchor="end">%1.1f</text>\n</g>', ...
                tc(i),tc(i),tc(i)+1,v(i));
        end
        v = linspace(timeline(1),timeline(end),6); tc = linspace(7.5,45,6);
        for i=1:6
            fh{end+1} = sprintf('<line x1="%1.4f" x2="%1.4f" y1="27" y2="26" stroke="black" stroke-width="0.2"/>\n<g transform="translate(%1.4f,29)">\n<text x="0" y="0" font-size="2" font-family="Verdana" text-anchor="middle">%1.1f</text>\n</g>', ...
                tc(i),tc(i),tc(i),v(i));
        end
    end
    fh{end+1} = sprintf('</g>\n<g transform="translate(72.5,68)">\n<text x="0" y="0" font-size="3" font-family="Verdana" text-anchor="middle">%s</text>\n</g>', ...
        ictitle);
else
    timec = axes('Position',[0.5,0.05,0.45,0.275]);
    image([timeline(1),timeline(end)],[mi,ma]+[1,-1]*(ma-mi)/4,max(min(cc,1),0),'Parent',timec);
    set(timec,'NextPlot','add');
    for i=1:size(timecourse,2)
        plot(timec,timeline,timecourse(:,i),'Color',[0;0;0],'LineWidth',2);
        set(timec,'NextPlot','add');
    end
    if (ma>mi), set(timec,'YLim',[mi,ma],'XLim',[timeline(1),timeline(end)]); end
    set(timec,'YDir','normal','XDir','normal','NextPlot','replace');
    set(get(timec,'Title'),'String',{ictitle});
end

if (outputsvg>1)
    fh = sprintf('<?xml version="1.0" encoding="UTF-8" standalone="no"?>\n<svg width="%1.4fpx" height="%1.4fpx" xmlns:svg="http://www.w3.org/2000/svg" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink">\n<g fill="none" stroke="black" stroke-width="0.5">\n<rect x="0" y="0" width="100" height="100" rx="2" ry="2"/>\n</g>%s</svg>', ...
        100,100,cat(2,fh{:}));
elseif (outputsvg>0)
    fh = sprintf('<g fill="none" stroke="black" stroke-width="0.5">\n<rect x="0" y="0" width="100" height="100" rx="2" ry="2"/>\n</g>%s', ...
        cat(2,fh{:}));
else
    set(fh,'Visible','on');
end

return


function data=scaleimg(data,newdim,f)

s = size(data);
if (length(s)==3)
    if (s(3)==3)
        nd = zeros([newdim,3]);
        for i=1:3
            nd(:,:,i) = scaleimg(squeeze(data(:,:,i)),newdim,f);
        end
        data = nd;
        return
    else
        return
    end
else
    [xx,yy] = meshgrid(linspace(1,size(data,2),newdim(2)), ...
        linspace(1,size(data,1),newdim(1)));
    if (f)
        data = interp2(data,xx,yy,'nearest');
    else
        data = interp2(data,xx,yy,'cubic');
    end
end


function a=mixit(a,b,f,o)

if (nargin<4), o = [1,1]; end
s = size(b);
f = repmat(f,[1,1,s(3)])*3/4;
a(o(1):(o(1)+s(1)-1),o(2):(o(2)+s(2)-1),:) = ...
    (a(o(1):(o(1)+s(1)-1),o(2):(o(2)+s(2)-1),:).*(1-f))+(b.*f);


function a=blowup(b,x,y)

s = size(b);
a = zeros(s.*[x,y]);
a(1:x:end,1:y:end) = b;


function a=rendererrorcomp(c,ec,sx,sy)

a = zeros(size(c).*[sx,sy]);
t = mfbox_mkgrid([sx,sy])-1;
l = sum((t-repmat(([sx,sy]-1)/2,size(t,1),1)).^2,2);
[l,s] = sort(l);
t = t(s,:);
ec = ceil(ec*sx*sy);
for i=1:size(t,1)
    a = a+circshift(blowup(c.*(ec>=i),sx,sy),[t(i,1),t(i,2)]);
end


function fh=makesvgimg(sx,sy,bb,c,ec,red,blue)

if (~isempty(bb))
    ttmp = tempname();
    bb = reshape(repmat(bb(:),1,3),[size(bb),3]);
    imwrite(max(min(bb,1),0),ttmp,'png');
    fh = sprintf('<image x="0" y="0" width="%d" height="%d" xlink:href="data:image/png;base64,%s"/>\n', ...
        sx,sy,mfbox_base64([],ttmp));
    delete(ttmp);
else
    fh = '';
end
t = zeros([numel(c),3]);
t(c>0,:) = repmat(red',sum(c(:)>0),1);
t(c<0,:) = repmat(blue',sum(c(:)<0),1);
sc = size(c);
p = c~=0 & ec~=0;
if (~isempty(bb))
    bb = scaleimg(bb,size(c),0);
    cc = reshape(mixit(bb,reshape(t,[size(c),3]),abs(c)),[],3);
else
    cc = t;
end
t = (mfbox_mkgrid(sc)-0.5)./repmat(sc,prod(sc),1);
t = t(p,:)*diag([sy,sx]);
cc = cc(p,:);
f = min(sx/sc(1),sy/sc(2));
ec = ec(p)*f;
cc = ceil(255*min(max(cc,0),1));
circ = cell(1,size(t,1));
for i=1:size(t,1)
    circ{i} = sprintf('<circle cx="%f" cy="%f" r="%f" fill="#%02x%02x%02x" stroke-width="0"/>', ...
        t(i,2),t(i,1),ec(i),cc(i,1),cc(i,2),cc(i,3));
end
fh = [fh,cat(2,circ{:})];


function bb=makepixelimg(bb,c,ec,red,blue)

c = rendererrorcomp(c,ec,4,4);
t = zeros([numel(c),3]);
t(c>0,:) = repmat(red',sum(c(:)>0),1);
t(c<=0,:) = repmat(blue',sum(c(:)<=0),1);
if (~isempty(bb))
    bb = scaleimg(bb,size(c),1);
    bb = reshape(repmat(bb(:),1,3),[size(bb),3]);
    bb = mixit(bb,reshape(t,[size(c),3]),abs(c));
else
    bb = mixit(ones([size(c),3]),reshape(t,[size(c),3]),abs(c));
end