TransmissionTomographyExample

Tiedoston sisältö

% Simulation code to study 
% transmission tomography problems (such as x-ray CT)
%
% Jyvaskyl Summer School 2014
% Aku Seppanen
% University of Eastern Finland, Kuopio
% Department of Applied Physics

clear all, close all
path(path,'.')

% Number of pixels in the true image: Nx1*Ny1;
Ny1 = 50;
%Ny1 = 120;
Nx1 = Ny1;

% Number of pixels in the reconstructed image: Nx*Ny;
Ny = 50;
Nx = Ny;

% Note: if Ny1 = Ny --> "Inverse crime"


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Create the true target image %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[XX1,YY1] = meshgrid(1:Nx1,1:Ny1);
XX1 = XX1(:);
YY1 = YY1(:);

% Target distribution with circular background and 2 inclusions
f_true = zeros(Nx1*Ny1,1);
cp1 = [Nx1/2,Nx1/2];
rad1 = [Nx1*2/5];
val1 = 1;
fnd1 = find((XX1-cp1(1)).^2 + (YY1-cp1(2)).^2 <= rad1^2);
f_true(fnd1) = val1;
cp2 = [Nx1/2,Nx1/4];
rad2 = [Nx1/15];
val2 = 3;
fnd1 = find((XX1-cp2(1)).^2 + (YY1-cp2(2)).^2 <= rad2^2);
f_true(fnd1) = val2;
cp3 = [Nx1/5,Nx1/2];
rad3 = [Nx1/10];
val3 = 2;
fnd3 = find((XX1-cp3(1)).^2 + (YY1-cp3(2)).^2 <= rad3^2);
f_true(fnd3) = val3;

F_true = reshape(f_true,Nx1,Ny1);
figure(1), clf, imagesc(F_true),axis image, colormap('hot'), colorbar
title('True image')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Simulate the observations(projection data corresponding to true target %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Minimum and maximum angle of projection
% note: in this geometry, if min_angle = 0, then max_angle < 180 means
% limited angle case.
min_angle = 0;
max_angle = 180; 
%max_angle = 90; 

% Number of projections
n_projections = 70;
%n_projections = 30;
%n_projections = 15;
%n_projections = 5;

% Projection angles
theta = linspace(min_angle,max_angle,n_projections);

% Simulate the projection data using Matlab function radon.m
[R] = radon(F_true,theta);
R = (1/Nx1)*R; % scaling by the length
% Check the number of pixels in the projection of an image in the inversion
% grid and interpolate the accurate projection data to this lower
% resolution pixel basis.
% (This is done to avoid inverse crimes.)
[R2] = radon(zeros(Nx,Ny),theta);
ng_theta = size(R2,1);
for ii = 1:size(R,2)
    R3(:,ii) = resample(R(:,ii),ng_theta,size(R,1));
end

% Noiseless data vector
g_true = R3(:);
ng = length(g_true);

% Variance of the observation noise: try different noise levels
var_n = .00001^2; % very small noise
%var_n = .001^2; % small noise
%var_n = .01^2; % not so small noise
%var_n = .1^2; % large noise
%Add noise to the data
n = sqrt(var_n)*randn(ng,1);
g = g_true + n;

figure(2),plot(g_true,'b'), hold on
plot(g,'r--')
legend('Noiseless observations','Noisy observations')

gg = reshape(g,ng/n_projections,n_projections);
figure(3), clf, imagesc([min_angle,max_angle],[0 ng_theta],gg,[min(g),max(g)]),axis image, colormap('hot'), colorbar
axis on, set(gca,'YDir','normal')
xlabel('Projection angle (degrees)')
ylabel('Sensor position in each projection')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Build up the observation matrix K %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[K,XX,YY] = ConstructObsMatrixTransTomo(Nx,Ny,theta);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Non-regularized LS solution %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f_LS = K\g;
F_LS = reshape(f_LS,Nx,Ny);
%figure, imagesc(F_LS),axis image, colormap('hot'), colorbar

minval = min([f_true;f_LS]);
maxval = max([f_true;f_LS]);
figure(4), clf
subplot(1,2,1), imagesc(F_true,[minval,maxval]),axis image, colormap('hot'), colorbar
title('True image')
subplot(1,2,2), imagesc(F_LS,[minval,maxval]),axis image, colormap('hot'), colorbar
title('LS solution')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Standard Tikhonov regularized solution %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

alpha = .01;
L_alpha = eye(Nx*Ny);
f_alpha = (K'*K + alpha*L_alpha'*L_alpha)\(K'*g);
F_alpha = reshape(f_alpha,Nx,Ny);
%figure, imagesc(F_alpha),axis image, colormap('hot'), colorbar

minval = min([f_true;f_alpha]);
maxval = max([f_true;f_alpha]);
figure(5), clf
subplot(1,2,1), imagesc(F_true,[minval,maxval]),axis image, colormap('hot'), colorbar
title('True image')
subplot(1,2,2), imagesc(F_alpha,[minval,maxval]),axis image, colormap('hot'), colorbar
title('Standard Tikhonov regularized solution')