Siv4 test.m code

From BanghamLab
Jump to navigation Jump to search

Return

This code illustrates lots of things to do with 1D sieving of 2D images column by column where the image is rotated.


function siv4_test
    %function siv4_test demonstrate siv4 functionality
    %close all
    figNo=0;
    subplotrows=[];
    if true
        
        %% Simple starting signal and a simplified signal
        figNo=newFig(figNo,'Simple starting signal and a simplified signal');
        data{1}=siv4_alt('PULSES3WIDE');
        showResult(subplotrows,...
            data{1}.name,data{1}.X,...
            sprintf('1D sieve columnwise scale %d',data{1}.options(1)),data{1}.y{1} ...
            );
        %% Showing increasing simplification by removing small scale pulses
        figNo=newFig(figNo,'Showing increasing simplification by removing small scale pulses');
        data{1}=siv4_alt('PULSES3WIDE',[2;5;10])
        showResult(4,...
            data{1}.name,data{1}.X,...
            sprintf('1D sieve columnwise scale %d',data{1}.options(1,1)),data{1}.y{1}, ...
            sprintf('1D sieve columnwise scale %d',data{1}.options(2,1)),data{1}.y{2}, ...
            sprintf('1D sieve columnwise scale %d',data{1}.options(3,1)),data{1}.y{3} ...
            );
        data{1}
%         figNo=newFig(figNo,'Showing increasing simplification by removing small scale pulses');
%         showResult(4,...
%             data{1}.name,data{1}.X,...
%             sprintf('1D sieve columnwise scale %d',data{1}.options(1,1)),data{1}.X-data{1}.y{1}, ...
%             sprintf('1D sieve columnwise scale %d',data{1}.options(2,1)),data{1}.X-data{1}.y{2}, ...
%             sprintf('1D sieve columnwise scale %d',data{1}.options(3,1)),data{1}.X-data{1}.y{3} ...
%             );

        %% Dynamic figure showing gradual simplification of signal by removing small scale details
        figNo=newFig(figNo,'Dynamic figure showing gradual simplification of signal by removing small scale details');
        im=getData('PULSES3WIDE');
        X(:,1)=double(im(:));
        maxmesh=size(X,1);
        data{1}=siv4_alt(X(:),maxmesh,'g');
        %     data{1}=siv4_alt('PULSES3WIDE',maxmesh,'g');
        cla
        shg
        plot(data{1}.X(:),'o-b');
        hold on
        g=data{1}.g{1};
        h=[];
        xlabel('position')
        ylabel('granule sample value')
        title('Granule values superimposed on data{1}');
        for i=2:length(g)
            if ~isempty(h) && ishandle(h)
                delete(h);
            end
            h=plot([g(1,i),g(1,i),g(1,i)+g(2,i),g(1,i)+g(2,i)],[0,g(3,i),g(3,i),0],'-r');
            pause(0.1)
        end
        
        %% Heatmap showing granules in scale-space
        figNo=newFig(figNo,'Heatmap showing granules in scale-space');
        cla
        maxHeat=max(g(4,2:end));
        minHeat=min(g(4,2:end));
        scales=unique(g(2,2:end));
        heatmap=zeros([length(X),max(scales)+2]);
        for i=1:length(scales)
            scale=scales(i);
            inds=find(g(2,:)==scale);
            for j=1:length(inds)
                start=g(1,inds(j));
                finish=start+g(2,inds(j))-1;
                heatmap(start:finish,scale)=heatmap(start:finish,scale)+g(4,inds(j))+3;
            end
        end
        heatmapR=(heatmap');
        imagesc(max(heatmapR(:))-heatmapR); % white background is nice
        xlabel('X')
        ylabel('Scale');
        colormap jet
        title('Heatmap, colormap jet')
        
        %% Aside: compare sieve filter bank (m-sieve) with Gaussian filter bank
        figNo=newFig(figNo,'Aside compare sieve filter bank (m-sieve) with Gaussian filter bank');
        clf
        scales=[0 1 2 4 8 16 32 64];
        for i=1:length(scales)
            scale=scales(i);
            if false %scale==0
                ys{i}=X;
                yg{i}=X;
            else
                ys{i}=siv4_alt(X,scale);
                h=fspecial('Gaussian',5*(scale+1),scale+0.1);
                midInd=round(5*(scale+1)/2);
                hh=h(midInd,:);
                hh=hh/sum(hh(:));
                yg{i}=conv(X',hh,'same');
            end
        end
        heatmapS=zeros([length(X),length(scales)]);
        heatmapG=zeros([length(X),length(scales)]);
        for i=1:length(scales)
            heatmapS(:,i)=cell2mat(ys{i}.y);
            heatmapG(:,i)=yg{i}';
        end
        subplot(1,2,1)
        imagesc(heatmapS')
        title('Sieve decomposition')
        str=;
        for i=1:length(scales)
            scale=scales(i);
            str=[str,num2str(scale),'|'];
        end
        set(gca,'Yticklabel',str)
        ylabel('Scale - note logarithmic');
        subplot(1,2,2)
        imagesc(heatmapG')
        title('Gaussian decomposition')
        colormap jet
        set(gca,'Yticklabel',str)
        ylabel('Scale - note logarithmic');
        
        %%  Histogram showing number of granules at each position (of various scales)
        figNo=newFig(figNo,'Histogram showing number of granules at each position (of various scales)');
        clf
        log_maxmesh=log(maxmesh); % want ten scale bins
        k=10/log_maxmesh;
        granularities=zeros(size(X,1),ceil(log_maxmesh*k));
        for i=2:length(g)
            for j=1:g(2,i)
                granularities(g(1,i)+j-1,1+ceil(log(g(2,i))*k))=granularities(g(1,i)+j-1,1+ceil(log(g(2,i))*k))+1;
            end
        end
        number_granules_at_each_position=sum(granularities,2);
        [maxnumber,i]=max(number_granules_at_each_position);
        bar(number_granules_at_each_position);
        title(sprintf('A stable extremum is at position %d with %d nested granules\n',i,maxnumber));
        xlabel('position (x)');
        ylabel('number of granules at each position');
        
        %%  Two dimensional signal with each column simplified to scale 10
        figNo=newFig(figNo,'Two dimensional signal with each column simplified to scale 10');
        data{1}=siv4_alt; % by default this reads in 'TestCard5.png' and sets all parameters
        showResult(subplotrows,...
            data{1}.name,data{1}.X,...
            sprintf('1D sieve columnwise scale %d',data{1}.options(1)),data{1}.y{1} ...
            );
        
        %%  Image simplified to scale 10 where columns are selected at angles
        figNo=newFig(figNo,'Image simplified to scale 10 column by column where the column is rotated');
        rotDegrees=[0 27 34 45]; % degrees from the vertical
        for i=1:length(rotDegrees)
            data{i}=data{1};
            if i==1
                [m,n]=size(data{1}.X);
                data{i}.scan=[m,m];
            else
                data{i}.scan=[5-i,1];
            end
            data{i}=siv4_alt(data{i}.X,data{i}.options,data{i}.outputs,data{i}.scan);
            subplot(2,2,i);
            imshow(uint8(data{i}.y{1}))
            title(sprintf('%d deg. scale %d, scan [%d %d]',rotDegrees(i),data{i}.options(1),data{i}.scan))
        end
        
        %%  Rotate: 1D sieve to scale 10 column by column: rotated back
        figNo=newFig(figNo,'Rotate: 1D sieve to scale 10 column by column: rotated back');
        rotDegrees=round(linspace(0,180,9)); % degrees from the vertical
        subplot(3,3,1);
        [~,requiredPlotParams]=EmbedAndImshow(data{1}.X,rotDegrees);
        subplot(3,3,1);
        data{9}.y{3}=EmbedAndImshow(uint8(data{1}.X),rotDegrees,requiredPlotParams,true);
        scaleToShow=data{1}.options(1);
        for i=1:length(rotDegrees)-1
            data{i}=data{1};
            rotDegree=rotDegrees(i);
            data{i}.X=imrotate(data{1}.X,rotDegree); % replace with rotated data
            data{i}=siv4_alt(data{i}.X,data{i}.options,data{i}.outputs,data{i}.scan);
            data{i}.y{2}=imrotate(data{i}.y{1},360-rotDegree,'crop');
            subplot(3,3,i+1);
            data{i}.y{3}=EmbedAndImshow(uint8(data{i}.y{2}),rotDegrees,requiredPlotParams,true);
            title(sprintf('%d deg. scale %d, scan [%d %d]',rotDegrees(i),data{i}.options(1),data{i}.scan))
        end
        
%         %%  Image granules from columns
%         figNo=newFig(figNo,'Image granules from columns');
%         data{1}=siv4_alt; % by default this reads in 'TestCard5.png' and sets all parameters
%         X=data{1}.X;
%         imSieved=siv4_alt(X,maxmesh,'g');
%         gtemp=cell2mat(imSieved.g);
%         g=gtemp(:,2:end); % omitting details on number of granules etc.
%         scales=1:8;
%         subplot(3,3,1)
%         imagesc(X)
%         title('Test image');
%         GRANULESasLINES=true;
%         for i=1:length(scales)
%             scale=scales(i);
%             ind=find(g(2,:)==scale);
%             granulesAtScale=g(:,ind);
%             subplot(3,3,i+1);
%             plotGranulesOnImage(X,granulesAtScale,scale,GRANULESasLINES);
%         end
%     else
        %%  Rotated image granules from columns rotated back
        verbose=true;
        viewRotated=true;
        
        figNo=newFig(figNo,'Rotated image granules from columns rotated back');
        data{1}=siv4_alt; % by default this reads in 'TestCard5.png' and sets all parameters
        %clf
        im=uint8(data{1}.X);
        rotDegrees=round(linspace(0,180,9)); % degrees from the vertical
        scalesToShow=[1,2,4,8,16,32,64];
        %scalesToShow=2;
        %GranFeaturesPerPixel=struct(
        for i=1:length(rotDegrees)-1
            data{i}=data{1};
            rotDegree=rotDegrees(i);
            % Rotate the image
            Rotim=imrotate(im,rotDegree,'loose'); % replace with rotated data
            [sievedRotim,gransRotim]=sieveRotated(Rotim,scalesToShow(end));
            for scale_i=2:length(scalesToShow)
                ind=find([gransRotim(2,1:end)]>=scalesToShow(scale_i-1) & ...
                    [gransRotim(2,1:end)]<scalesToShow(scale_i));
                if ~isempty(ind)
                    temp_gransRotim=[gransRotim(:,1), gransRotim(:,ind)];
                    scaleToShow=scalesToShow(scale_i);
%                     figure(1); colormap gray
%                     cla; hold off
%                     imagesc(double(Rotim)); hold on;
                    title(sprintf('%d',rotDegree));
                    % Rotate the sieved image back
                    figure(figNo); colormap gray
                    cla; hold off
                    [GransRotBack]=rotate_sieved_image_back(...
                        im,size(Rotim),rotDegree,temp_gransRotim,...
                        scaleToShow,verbose,viewRotated,Rotim);
                    disp('press a key')
                    pause(0.2)
                end
            end
        end
    end
end
%%
function [sievedRotim,gransRotim]=sieveRotated(Rotim,scaleToShow)
    filename=sprintf('sievedData-%s.mat',num2str(scaleToShow));
%     if exist(filename)==2
%         load(filename)
%     else
        data=siv4_alt(double(Rotim),scaleToShow,'l'); % rotated and lowpass sieved
        sievedRotim=data.y{1};
        SIVedGData=siv4_alt(double(Rotim),scaleToShow,'g'); % rotated and sieved to granules
        gransRotim=cell2mat(SIVedGData.g);
        save(filename, 'sievedRotim','gransRotim')
%     end
end
%%

function [GransRotBack]=rotate_sieved_image_back(...
        im,sizRotim,rotDegree,gransRotim,gran_scale,verbose,viewRotated,Rotim)
    % [RotImBack,GransRotBack]=rotate_sieved_image_back(im,imIn,sizRotim,rotDegree,gransRotim,gran_scale,padding,verbose)
    %
    % RotImBack,  Rotated imIn (back) by 360-rotDegrees
    % GransRotBack, granules computed from the rotated image 
    %          (gransRotim) are also rotated back (i.e. the first row - 
    %           ignoring the first column which means something different)
    if nargin<8
        verbose=true;
    end
    if nargin<9
        viewRotated=true;
    end
    %RotImBack=imrotate(imIn,360-rotDegree,padding);

    Yin=sizRotim(1);
    Xin=sizRotim(2);
    [Yout,Xout]=size(im); % rotated size
    %[Yout,Xout]=size(RotImBack); % rotated size
    if verbose
        if viewRotated
            imagesc(double(Rotim)); hold on;
        else
            imagesc(double(im)); hold on;
        end
    end
    GransRotBack.gransRotim=gransRotim;
    GransRotBack.gransRotim(1,2:end)=0; % clear the indices - not all will be filled.
    gran_scale=min(gran_scale,max(gransRotim(2,2:end)));
    % Rotate the granules back
    for j=2:size(gransRotim,2)
        granRotim=gransRotim(:,j);
        [yyr,xxr]=ind2sub([Yin,Xin],granRotim(1));
        yyr_start=yyr+0.5*gran_scale;
        yyr_end=yyr-0.5*gran_scale;
        % now rotate these back, start point
        backRotDegree=360-rotDegree;
        xyr=ceil([cosd(backRotDegree),sind(backRotDegree);-sind(backRotDegree),...
            cosd(backRotDegree)]*[xxr-Xin/2;yyr-Yin/2]+[Xout/2;Yout/2]);
        %xx=min(max(xyr(1),1),Xin);
        %yy=max(min(xyr(2),Yin),1);
        xx=min(max(xyr(1),1),Xout);
        yy=max(min(xyr(2),Yout),1);
        GransRotBack.gransRotim(1,j)=sub2ind([Yout,Xout],yy,xx);
        % and the end point
        yy_start=yy+0.5*gran_scale;
        yy_start=max(min(yy_start,Yout),1);
        yy_end=yy-0.5*gran_scale;
        yy_end=max(min(yy_end,Yout),1);
        yy_range=ceil(yy_end:yy_start);
        GransRotBack.indexList{j}=sub2ind([Yout,Xout],yy_range,xx*ones(size(yy_range)));
        if verbose
            if viewRotated
                title(sprintf('scale %d',gran_scale))
                xyer=[xyr(1);xyr(2)-1];
                plot([xxr xxr],[yyr_start yyr_end],'-r');
                plot(xxr,yyr,'*r');
            else
                xyer=[xyr(1);xyr(2)-1];
                plot([xx xx],[yy_start yy_end],'-r');
                plot(xx,yy,'*r');
            end
        end
    end
end

%%
function [X,Y]=EmbedAndPlot(xy,im,rotDegrees,requiredPlotParams,verbose)
    if nargin<5
        verbose=true;
    end
    X=xy(1);Y=xy(2);
    [rows,cols]=size(im);
    offsetI=floor((requiredPlotParams.rowsMax-rows)/2);
    offsetJ=floor((requiredPlotParams.colsMax-cols)/2);
    if size(xy,2)>1
        X=xy(1,:)+offsetJ+1;
        Y=xy(2,:)+offsetI-1;
        if verbose
            hold on
            y=repmat(requiredPlotParams.rowsMax,1,size(Y,2))-Y;
            plot(X,y,'-b');
        end
    else
        X=xy(1)+offsetJ+1;
        Y=xy(2)+offsetI-1;
        if verbose
            hold on
            plot(X,requiredPlotParams.rowsMax-Y,'.b');
        end
    end
    
end

%%
function plotGranulesOnImage(X,granulesAtScale,scale,GRANULESasLINES)
    %subplot(1,2,2)
    imagesc(X)
    title(sprintf('Scale %d column granules',scale));
    hold on
    colormap gray
    im=zeros(size(X));
    for i2=1:length(granulesAtScale)
        G=granulesAtScale(:,i2);
        [ii,jj]=ind2sub(size(X),G(1));
        if GRANULESasLINES
            if G(4)>0
                if scale==1
                    plot(jj,ii,'*r','linewidth',2);
                else
                    plot([jj, jj],[ii ii+scale],'-r','linewidth',2);
                end
            else
                if scale==1
                    plot(jj,ii,'*b','linewidth',2);
                else
                    plot([jj, jj],[ii ii+scale],'-b','linewidth',2);
                end
            end
        else
            for k=1:scale
                if G(4)>0
                    plot(jj,ii+k-1,'.r');
                else
                    plot(jj,ii+k-1,'.b');
                end
                im(ii:ii+k,jj)=G(4);
            end
        end
    end
    pause(1)
end

function figNo=newFig(figNo,titleString)
    screensize=get(0,'screensize');
    figNo=figNo+1;
    figure(figNo);
    pos=get(figNo,'Position');
    pos(1)=10+figNo*15;
    if figNo==1
        pos(2)=screensize(4)-pos(4)-90*figNo;
    else
        pos(2)=screensize(4)-pos(4)-60*figNo;
    end
    set(figNo,'position',pos,'name',titleString);
end


%%
function [tempX,requiredPlotParams]=EmbedAndImshow(X,rotDegrees,requiredPlotParams,verbose)
    if nargin<4
        verbose=true;
    end
    tempX=X;
    if nargin<3 % we are creating params from X
        maxEdgeLength=max(size(X));
        for i=1:length(rotDegrees)-1
            rotDegree=rotDegrees(i);
            XR=imrotate(X,rotDegree); % imshow will then be consistent with padding
            sXR=size(XR);
            if max(sXR(:))>maxEdgeLength
                maxEdgeLength=max(sXR);
            end
        end
        requiredPlotParams.ax=[0,maxEdgeLength,0,maxEdgeLength];
        [requiredPlotParams.rowsX,requiredPlotParams.colsX]=size(X);
        requiredPlotParams.rowsMax=maxEdgeLength;
        requiredPlotParams.colsMax=maxEdgeLength;
    else% we are using requiredPlotParams to imshow the image X
        tempX=zeros([requiredPlotParams.rowsMax,requiredPlotParams.colsMax]);
        [rows,cols]=size(X);
        offsetI=floor((requiredPlotParams.rowsMax-rows)/2);
        offsetJ=floor((requiredPlotParams.colsMax-cols)/2);
        tempX(1+offsetI:offsetI+rows,1+offsetJ:offsetJ+cols)=X;
        %         tempX=zeros([requiredPlotParams.rowsMax,requiredPlotParams.colsMax]);
        %         for j=1:cols %requiredPlotParams.cols45-offsetJ-1
        %             tempX((offsetI+1:requiredPlotParams.rowsMax-offsetI)',j+offsetJ)=X(:,j);
        %         end
        if verbose
            title(sprintf('rotated by %d degrees',rotDegrees));
            imshow(uint8(tempX));
        end
    end
end

%%
function showResult(vcols,varargin)
    if ~isempty(vcols)
        vrows=1;
        k=1;
        for j=1:vcols
            subplot(vrows,vcols,k)
            result=varargin{k*2};
            if k==1
                signal=result;
            end
            if any(size(result)==1) % then it is 1D
                plot(result,'o-b');
                hold on
                plot(signal,':b');
                hold off
            else
                imshow(uint8(result))
            end
            title(varargin{k*2-1});
            k=k+1;
        end
    else
        v_elements=numel(varargin)/2;
        vrows=floor(sqrt(v_elements));
        vcols=ceil(v_elements/vrows);
        k=1;
        for i=1:vrows
            for j=1:vcols
                subplot(vrows,vcols,k)
                result=varargin{k*2};
                if k==1
                    signal=result;
                end
                if any(size(result)==1) % then it is 1D
                    plot(result,'o-b');
                    hold on
                    plot(signal,':b');
                    hold off
                else
                    imshow(uint8(result))
                end
                title(varargin{k*2-1});
                k=k+1;
            end
        end
    end
end