Siv4 test.m code
Jump to navigation
Jump to search
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