clc;

clear;

close all;

warning off;

pack;

addpath 'func\'

%Pixel Size

Pix_Size = 32;

[I,E] = phantom(Pix_Size);

figure;

imshow(I);

[rays,sino] = siddon(I);

ind = find(sum(rays,2));

A = rays(ind,:);

b = sino(:);

b = b(ind);

%calculate the singular value of A

s = svds(A,size(A,2));

%plot loglog figure

figure;

subplot(121);

plot(s,'b-o');

axis([0,size(A,2),0,60]);

grid on;

axis square;

subplot(122);

loglog(s,'b-o');

axis([0,size(A,2),0,60]);

grid on;

axis square;

%Delete rows from the matrix and comment on the effect of this on the singular values

%Delete 500 rows

A2 = A;

A2(1:500,:) = [];

s2 = svds(A2,size(A2,2));

%Delete 1000 rows

A3 = A;

A3(1:1000,:)= [];

s3 = svds(A3,size(A3,2));

%Delete 2000 rows

A4 = A;

A4(1:2000,:) = [];

s4 = svds(A4,size(A4,2));

figure;

plot(s,'b');

hold on;

plot(s2,'r');

hold on;

plot(s3,'k');

hold on;

plot(s4,'g');

hold off;

axis([0,size(A,2),0,60]);

grid on;

legend('Initial singular values','singular values after delete 500 rows'

clc;

clear;

close all;

warning off;

pack;

addpath 'func\'

%Pixel Size

Pix_Size = 32;

I = phantom(Pix_Size);

figure;

subplot(121);

imshow(I);

title('Images');

%calculate A and b

[rays,sino] = siddon(I);

ind = find(sum(rays,2));

A = rays(ind,:);

b = sino(:);

b = b(ind);

%backslash operator means when AX = B then X=A/B

I2 = A\b;

%get the image matrix

I22 = reshape(I2,32,32);

subplot(122);

imshow(I22);

title('Reconstruct Images');

PSNR = psnr(I,I22)

%add error to b with 0.01 power

b2 = b + 0.01*randn(size(b));

I2 = A\b2;

I23 = reshape(I2,32,32);

%calculate the psnr

PSNR = psnr(I,I23)

%add error to b with 0.1 power

b2 = b + 0.1*randn(size(b));

I2 = A\b2;

I24 = reshape(I2,32,32);

%calculate the psnr

PSNR = psnr(I,I24)

%add error to b with 1 power

b2 = b + randn(size(b));

I2 = A\b2;

I25 = reshape(I2,32,32);

%calculate the psnr

PSNR = psnr(I,I25)

figure;

subplot(131);

imshow(I23);

title('err = 0.01');

subplot(132);

imshow(I24);

title('err = 0.1');

subplot(133);

imshow(I25);

title('err = 1');

%the function of calculating PSNR

function PSNR = psnr(f1, f2)

k = 8;

fmax = 2.^k - 1;

a = fmax.^2;

e = double(f1) - double(f2);

[m, n] = size(e);

b = sum(sum(e.^2));

PSNR = 10*log(m*n*a/b);