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