Siv4 test.m code

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