Radon transform applied to lines detection

Contents

DATA GENERATION

W=65;               % image width (must be odd W=2M+1)
H=65;               % image height
K=3;                % number of lines
spread=1;           % std of the gaussian blur
noiselevel=0;       % level of noise
randomgen=0;        % boolean if lines are generated randomly of manually
rng(0);             % seed of random numbers generator
plotData=1;         % display images v, x and y
plotComp=0;         % display others comparaisons theo vs. emp

t_k=[pi/6 -pi/5 pi/16]; % array containing angles of lines
a_k=[255 255 255];      % array containing amplitude of lines
p_k=[10 0 -15];         % array containing offset of lines

data_generation;    % Generate the blurred image xstar
                    % of these K lines with additional
                    % noise y=xstar+randn(H,W)*noiselevel

RADON TRANSFORM OF XSTAR (NORMAL BLUR)

Radon transform on the noisy data y

alpha = -90+0.414:0.1:90; % +0.414 to avoid the samples to correspond
[R,xp] = radon(y,alpha);  % to exact degrees of the lines t_k

Decomment to find the peaks (automaticaly) via function FastPeakFind

% Find peaks of the Radon transform (automaticaly)
% p=FastPeakFind(R);
% we remove this false peak
% p=[p(1:2);p(5:end)];
% palpha=p(1:2:end);
% pxp=p(2:2:end);

Find peaks of the Radon transform (manually)

pxp=[48 57 33];         % Vertical coordinates of the peaks
[M1 I1]=max(R(48,:));
[M1 I2]=max(R(57,:));
[M1 I3]=max(R(33,:));
palpha=[I1 I2 I3];      % Horizontal coordinates of the peaks

Peaks of the Radon transform - Estimation of lines parameters

alphapeaks=alpha(palpha)*pi/180;
xppeaks=xp(pxp)';
[t_k_emp,jj]=sort(alphapeaks);
[t_k,ii]=sort(t_k);
p_k_emp=xppeaks./cos(alphapeaks);

Display the estimated parameters and the relative errors

disp(['t_k = ',num2str(t_k)]);
disp(['t_k_emp = ',num2str(t_k_emp)]);
disp(['p_k = ',num2str(p_k)]);
disp(['p_k_emp = ',num2str(p_k_emp)]);
disp(['(t_k-t_k_emp)/t_k = ',num2str((t_k-t_k_emp)./t_k)]);
disp(['p_k-p_k_emp = ',num2str(p_k(ii)-p_k_emp(jj))]);
t_k = -0.62832     0.19635      0.5236
t_k_emp = -0.62982     0.19398     0.52035
p_k = 10   0 -15
p_k_emp = 0      10.3729     -15.2867
(t_k-t_k_emp)/t_k = -0.0023889    0.012089      0.0062
p_k-p_k_emp = 0     0.28669    -0.37292

Plot the Radon transform + peaks

figure, imagesc(alpha, xp, R); colormap(hot);
xlabel('$$\alpha$$ (degrees)','Interpreter','latex');
ylabel('$$x\prime$$','Interpreter','latex');
title('$$R_{\alpha} (x\prime)$$','Interpreter','latex');
colorbar
set(gcf,'color','w');
axis square;
hold on;
plot(alpha(palpha),xp(pxp),'gs')

Study the variation of relative error with respect to alpha discretization of the Radon transform

Tabstep=[1 0.1 0.01 0.001 0.0001];
NbSteps=length(Tabstep);
TabErr1=zeros(1,NbSteps);
TabErr2=zeros(1,NbSteps);
TabErr3=zeros(1,NbSteps);
for i=1:NbSteps
alpha = -90+0.414:Tabstep(i):90;
[R,xp] = radon(y,alpha);
[M I1] = max(R(48,:));
[M I2] = max(R(57,:));
[M I3] = max(R(33,:));
thetatilde1=alpha(I1)*pi/180;
thetatilde2=alpha(I2)*pi/180;
thetatilde3=alpha(I3)*pi/180;
TabErr1(i)=abs((t_k(1)-thetatilde1)/t_k(1));
TabErr2(i)=abs((t_k(3)-thetatilde2)/t_k(3));
TabErr3(i)=abs((t_k(2)-thetatilde3)/t_k(2));
end
figure;
semilogx(Tabstep,TabErr1,'d-b','MarkerFaceColor','b','MarkerSize',10);
hold on;
semilogx(Tabstep,TabErr2,'s-r','MarkerFaceColor','r','MarkerSize',10);
semilogx(Tabstep,TabErr3,'^-m','MarkerFaceColor','m','MarkerSize',10);
xlabel('Angle step');
ylabel('Relative Error');
title('Error estimation of each line''s angle respect to radon''step angle');
set(gcf,'color','w');
legend({'$$\theta=-\frac{\pi}{5}$$','$$\theta=\frac{\pi}{6}$$',...
    '$$\theta=\frac{\pi}{16}$$'},'Interpreter','latex','Location','NorthWest');

RADON TRANSFORM OF Y (STRONG NOISE)

noiselevel=200;       % level of noise
t_k=[pi/6 -pi/5 pi/16]; % array containing angles of lines
a_k=[255 255 255];      % array containing amplitude of lines
p_k=[10 0 -15];         % array containing offset of lines

data_generation;    % Generate the blurred image xstar
                    % of these K lines with additional
                    % noise y=xstar+randn(H,W)*noiselevel

Radon transform on the noisy data y

alpha = -90+0.414:0.1:90; % +0.414 to avoid the samples to correspond
[R,xp] = radon(y,alpha);  % to exact degrees of the lines t_k

Plot the Radon transform + peaks

figure, imagesc(alpha, xp, R); colormap(hot);
xlabel('$$\alpha$$ (degrees)','Interpreter','latex');
ylabel('$$x\prime$$','Interpreter','latex');
title('$$R_{\alpha} (x\prime)$$','Interpreter','latex');
colorbar
set(gcf,'color','w');
axis square;
hold on;
plot(alpha(palpha),xp(pxp),'gs')

RADON TRANSFORM OF XSTAR (STRONG BLUR)

spread=8;           % std of the gaussian blur
noiselevel=0;       % level of noise
t_k=[pi/6 -pi/5 pi/16]; % array containing angles of lines
a_k=[255 255 255];      % array containing amplitude of lines
p_k=[10 0 -15];         % array containing offset of lines

data_generation;    % Generate the blurred image xstar
                    % of these K lines with additional
                    % noise y=xstar+randn(H,W)*noiselevel

Radon transform on the noisy data y

alpha = -90+0.414:0.1:90; % +0.414 to avoid the samples to correspond
[R,xp] = radon(y,alpha);  % to exact degrees of the lines t_k

Plot the Radon transform + peaks

figure, imagesc(alpha, xp, R); colormap(hot);
xlabel('$$\alpha$$ (degrees)','Interpreter','latex');
ylabel('$$x\prime$$','Interpreter','latex');
title('$$R_{\alpha} (x\prime)$$','Interpreter','latex');
colorbar
set(gcf,'color','w');
axis square;
hold on;
plot(alpha(palpha),xp(pxp),'gs')