CT图像重建算法——联合代数重建算法(SART)

联合代数重建算法(simultaneous ART,SART)是A.H.Andersen和A.C.Kak于1984年提出一种改进的迭代重建算法，只需要很少的选代次数就可以得到很好的重建质量和精度。
ART算法每次迭代只用到一条射线投影,测量噪声很容易被引人，并且需要较大的迭代次数，重建效率不高。SART算法对ART算法做了一些改进：利用同一投影角度下通过像素的所有射线的误差来确定对该个像素进行校正值，而不是只考虑一条射线。这样就相当于对ART中的噪声进行了平滑，重建图像对测量噪声就不那么敏感了。SART的迭代公式为

可以看到，SART不是按照逐条射线对图像像素进行更新，而是在计算完一个特定投影角度的整个投影之后再进行更新。
SART算法与另外一种经典算法SIRT(Simultaneous lterative Reconstruction Technique .联合迭代重建算法)非常相似，二者都属于联合迭代形式，都能抑制测量误差和一些干扰因素。不同的是，在一次迭代更新中SART算法使用的是某一投影角度下通过某一像素的所有射线.而SIRT算法用到的是通过该像素的所有射线。SIRT算法有一个很大的缺点--收敛速度慢以及重建时间很长，这导致其并没有得到广泛的普及。SART被认为是结合了ART和SIRT两种算法所有的优点而丢弃了它们的缺点。

clc;
clear all;
close all;
N=180;%原始
N2=N^2;
I=phantom(N);
theta=linspace(0,180,61);
theta=theta(1:60);
%产生投影数据
P_num=260;%探测器通道
P=ParallelBFP(theta,N,P_num);%解析法产生投影数据
% P=radon(I,theta);
%获取投影矩阵
delta=1;
[W_ind,W_dat]=SystemMatrixV2(theta,N,P_num,delta);
%设置参数
F=zeros(N2,1);
lamda=0.25;%松弛因子 
c=0;
irt_num=5;%总迭代次数 
while (c<irt_num)
    for j=1:length(theta)
        W1_ind=W_ind((j-1)*P_num+1:j*P_num,:);
        W1_dat=W_dat((j-1)*P_num+1:j*P_num,:);
        W=zeros(P_num,N2);
        for jj=1:P_num
            %如果射线不经过任何像素，不作计算
            if ~any(W1_ind(jj,:))
                continue
            end
            for ii=1:2*N
                m=W1_ind(jj,ii);
                if m>0 && m<=N2
                    W(jj,m)=W1_dat(jj,ii);
                end
            end
        end
        sumCol=sum(W)';%列和向量
        sumRow=sum(W,2);%行和向量
        ind1=find(sumRow>0);
        corr=zeros(P_num,1);
        err=P(:,j)-W*F;
        corr(ind1)=err(ind1)./sumRow(ind1);%修正误差
        backproj=W'*corr;%修正误差反投影
        ind2=find(sumCol>0);
        delta=zeros(N2,1);
        delta(ind2)=back