Radon变换的matlab仿真

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','singular values after delete 1000 rows','singular values after delete 2000 rows');

 

fig1. 32*32 Shepp-Logan phantom

 

fig2. singular value common plot(left) and log-log plot(right)

 

fig3. singular value after deleting some 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);

 

fig4. the initial image and the reconstruct images

     The PSNR is 795,which mean the quality of reconstruct image is very good.

fig5. different quality of different noise power