matlab_医学CT重建 ART，SART算法

clear all close all tic N = 256; n = 0;%迭代次数 it0 = ones(4*N);%数据矩阵扩增为原来四倍，保证对旋转过程的灵敏度 it1=it0; sl=(0:1:100); slf=(0:1:99); angle=(0:180)'; nangle = length(angle); I = phantom('Modified Shepp-Logan',256); %I=imread('C:\Users\adminraden\Desktop\医学成像系统大作业\CT图像\800px-CT_ScoutView.jpg'); %I=rgb2gray(I); %I=double(I); %I = imresize(I,[256 256]); %I=I./255; P = radon (I, angle); R0 = iradon (P, angle); R0 = imresize(R0,[256 256]); P1 = radon (I, 0:10:180); %P1=P; %P1(200,1:181)=100; %P1(200,1:181)=(P1(199,1:181)+P1(201,1:181))/2; R1 = iradon (P1, 0:10:180); %R1 = iradon (P1, angle); for a = 1:N for b = 1:N it0 ((4*a-3):4*a, (4*b-3):4*b) = R1 (a, b)/16; I1 ((4*a-3):4*a, (4*b-3):4*b) = I (a, b)/16; end end while (n < 100) n = n + 1; for ii = 1:N pj = zeros (4*N); pj (:,(4*ii-3):4*ii) = 1; pj1 = imrotate (pj, angle (mod(n-1, nangle)+1), 'crop');%旋转采样 P_pj = sum(sum(I1.*pj1)); R_pj = sum(sum(it0.*pj1)); it1 = it0 + pj1.*(P_pj-R_pj)/sum(sum(pj1)); it0 = it1; end % for ii = N:2 % pj = zeros (4*N); % pj (:,(4*ii-3):4*ii) = 1; % pj1 = imrotate (pj, angle (mod(n-1, nangle)+1), 'crop');%旋转采样 % % P_pj = sum(sum(I1.*pj1)); % R_pj = sum(sum(it0.*pj1)); % % it1 = it0 + pj1.*(P_pj-R_pj)/sum(sum(pj1)); % it0 = it1; % end for a = 1:N for b = 1:N R (a, b) = sum (sum (it0 ( (4*a-3):4*a, (4*b-3):4*b))); end end dA=R-R0; % slf(n)=max(max(abs(dA))); slf(n)=norm(dA,1); end % for a = 1:N; % for b = 1:N; % R (a, b) = sum (sum (it0 ( (4*a-3):4*a, (4*b-3):4*b))); % end % end figure, subplot(2,2,1);imshow(I), title('原始图像');stem(slf); subplot(2,2,2);imshow(R0), title('完整的采集信号重构图像'); subplot(2,2,3);imshow(R1), title('不完整信号重构图像'); subplot(2,2,4);imshow(R), title('迭代重建图像'); toc