% 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')