ART算法

[兮追儿]

clc;
clear all;
close all;
N = 180;
N2 = N^2;
I = phantom(N);
theta = linspace(0,180,181);
theta = theta(1:180);
P_num = 260;
P = radon(I,theta);
delta = 1;
[W_ind, W_dat] = medfuncSystemMatrix(theta, N, P_num, delta);
F = zeros(N2,1);
lambda = 0.25;
c = 0;
irt_num = 10;
while(c<irt_num)
for j=1:length(theta)
for i=1:1:P_num
u = W_ind((j-1)*P_num+i,:);
v = W_dat((j-1)*P_num+i,:);
if any(u) == 0
continue;
end
w = zeros(1,N2);
ind = u>0;
w(u(ind)) = v(ind);
PP = w*F;
C = (P(i,j) - PP)/sum(w.^2)*w';
F = F + lambda*C;
end
end
F(F<0) = 0;
c = c + 1;
end
F = reshape(F,N,N)';
figure, hold on
subplot(1,2,1),imshow(I),xlabel('(a)256*256 head phantom');
subplot(1,2,2),imshow(F),xlabel('(b)ART reconstructed head phantom');