CT图像重建算法——迭代算法ART

迭代图像重建的方法可分为迭代迭代法和统计迭代法，代数迭代法以代数重建算法(Algebraic Reconstruction Technique,ART)为代表。ART适合于不完全投影数据的图像重建，其抗噪声干扰能力强，另外可以结合一些先验知识进行求解，ART最大的缺点是计算量大，重建速度慢。

代数重建算法(ART)是由Kaczmarz于1937年在求解相容线性方程组时提出的，随后得到 Tanabe 的进一步阐明。ART先假设一幅初始图像，通过“正向投影”得到一个计算的投影，然后将计算的投影与实际测量投影相比较，基于比较的差值来估计修正值，将修正值均匀地分配给射线经过的那些像素上，逐条射线地执行这一过程，直到满足要求。

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);
%进行ART迭代
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
          w=zeros(1,N2);
          ind=find(u>0);
          w(u(ind))=v(ind);
          PP=w*F;
          error=P(i,j)-PP;
          C=error/sum(w.^2).*w';
          F=F+lambda*C;
      end
    end
    F(F<0)=0;%小于0的像素值置为0；
    c=c+1;
end
F=reshape(F,N,N)';
figure(1);
imshow(I);xlabel('(a)180*180头模型图像');
figure(2);
imshow(F);xlabel('(b)ART算法重建的图像');





%计算投影矩阵
function [W_ind,W_dat]=SystemMatrixV2(theta,N,P_num,delta)
%P_num:每个投影角度下的射线条数（探测器通道个数)
%delta:网格大小
%W_ind:存储投影射线所穿过的网格的编号的矩阵，M行，2*N列
%W_dat:存储投影射线所穿过的网格的长度的矩阵，M行，2*N列

%%用于验证的一小段程序
% theta=45;
% N=10;
% P_num=15;
% delta=1;
%------
N2=N^2;
M=length(theta)*P_num;%投影射线的总条数
W_ind=zeros(M,2*N);
W_dat=zeros(M,2*N);
% t_max=sqrt(2)*N*delta;
% t=linspace(-t_max/2,t_max/2,P_num);
t=(-(P_num-1)/2:(P_num-1)/2+1);
%t=(-(P_num-1)/2:(P_num-1)/2)*delta;%探测器坐标

for jj=1:length(theta)
    th=theta(jj);
    for ii=1:1:P_num
        %完成一条射线权因子向量的计算
        u=zeros(1,2*N);
        v=zeros(1,2*N);
       
        if th==0
            %如果超出网格范围，计算下一条
            if (t(ii)>=N/2*delta) || (t(ii)<=-N/2*delta)
                continue
            end
            kin=ceil(N/2+t(ii)/delta);
            kk=kin:N:(kin+N*(N-1));
            u(1:N)=kk;
            v(1:N)=ones(1,N)*delta;
           
            
        elseif th==90
            %如果超出网格范围，计算下一条
            if (t(ii)>=N/2*delta) || (t(ii)<=-N/2*delta)
                continue
            end
            kout=N*ceil(N/2-t(ii)/delta);
            kk=(kout-(N-1)):kout;
            u(1:N)=kk;
            v(1:N)=ones(1,N)*delta;
            %%投影角度等于0时%%
       
        else
            if th>90
                th_temp=th-90;
            elseif th<90
                th_temp=90-th;
            end
            th_temp=th_temp*pi/180;
            %确定射线y=mx+b的m和b
            b=t/cos(th_temp);
            m=tan(th_temp);
            y1d=(-N/2)*delta*m+b(ii);
            y2d=(N/2)*delta*m+b(ii);
            if (y1d<-N/2*delta && y2d<-N/2*delta)||(y1d>N/2*delta && y2d>N/2*delta)
                continue
            end
            %%确定入射点（xin,yin)、出射点（xout,yout)及参数d1%
            if y1d<=N/2*delta && y1d>=-N/2*delta && y2d>N/2*delta
                yin=y1d;
                yout=N/2*delta;
                xout=(yout-b(ii))/m;
                kin=N*floor(N/2-yin/delta)+1;
                kout=ceil(xout/delta)+N/2;
                d1=yin-floor(yin/delta)*delta;             
            elseif y1d<=N/2*delta && y1d>=-N/2*delta && y2d>=-N/2*delta && y2d<N/2*delta
                yin=y1d;
                yout=y2d;
                kin=N*floor(N/2-yin/delta)+1;
                kout=N*floor(N/2-yout/delta)+N;
                d1=yin-floor(yin/delta)*delta;
            elseif y1d<-N/2*delta && y2d>N/2*delta
                yin=-N/2*delta;
                yout=N/2*delta;
                xin=(yin-b(ii))/m;
                xout=(yout-b(ii))/m;
                kin=N*(N-1)+N/2+ceil(xin/delta);
                kout=ceil(xout/delta)+N/2;
                d1=N/2*delta+(floor(xin/delta)*delta*m+b(ii));
            elseif y1d<-N/2*delta && y2d>=-N/2*delta && y2d<N/2*delta
                yin=-N/2*delta;
                yout=y2d;
                xin=(yin-b(ii))/m;
                kin=N*(N-1)+N/2+ceil(xin/delta);
                kout=N*floor(N/2-yout/delta)+N;
                d1=N/2*delta+(floor(xin/delta)*delta*m+b(ii));
            else
                continue
            end
            %计算射线i穿过像素的编号和长度%
            k=kin;
            c=0;
            d2=d1+m*delta;
            while k>=1 && k<=N2
                c=c+1;
                if d1 >=0 && d2>delta
                    u(c)=k;
                    v(c)=(delta-d1)*sqrt(m^2+1)/m;
                    if k>N && k~=kout
                        k=k-N;
                        d1=d1-delta;
                        d2=d1+m*delta;
                    else
                        break;
                    end
                elseif d1>=0 && d2==delta
                    u(c)=k;
                    v(c)=delta*sqrt(m^2+1);
                    if  k>N && k~=kout
                        k=k-N-1;
                        d1=0;
                        d2=d1+m*delta;
                    else
                        break
                    end
                elseif d1>=0 && d2<delta
                    u(c)=k;
                    v(c)=delta*sqrt(m^2+1);
                    if k~=kout
                        k=k+1;
                        d1=d2;
                        d2=d1+m*delta;
                    else
                        break
                    end
                elseif d1<=0 && d2>=0 && d2<=delta
                    u(c)=k;
                    v(c)=d2*sqrt(m^2+1)/m;
                    if k~=kout
                        k=k+1;
                        d1=d2;
                        d2=d1+m*delta;
                    else
                        break
                    end
                elseif d1<=0 && d2>delta
                    u(c)=k;
                    v(c)=delta*sqrt(m^2+1)/m;
                    if k>N && k~=kout
                        k=k-N;
                        d1=d1-delta;
                        d2=d1+m*delta;
                    else
                        break
                    end
                elseif d1<=0 && d2==delta
                    u(c)=k;
                    v(c)=delta*sqrt(m^2+1)/m;
                    if k>N && k~=kout
                        k=k-N-1;
                        d1=0;
                        d2=m*delta;
                    else
                        break             
                    end
                end
              
            end
            %如果投影角度小于90，还需要利用投影射线关于y轴的对称性计算出权因子向量
            if th<90 
                u_temp=zeros(1,2*N);
                if any(u)==0
                    continue;
                end
                ind=find(u>0);
                for k=1:length(u(ind))
                    r=rem(u(k),N);
                    if r==0
                        u_temp(k)=u(k)-N+1;
                    else
                        u_temp(k)=u(k)-2*r+N+1;
                    end
                end
                u=u_temp;
            end
        end
        W_ind((jj-1)*P_num+ii,:)=u;
        W_dat((jj-1)*P_num+ii,:)=v;
    end
end

function P=ParalleBFP(theta,N,N_d)
shep=[1 .69 .92 0 0 0
    -.8 .6624 .8740 0 -.0184 0
    -.2 .1100 .3100 .22 0 -18
    -.2 .1600 .4100 -.22 0 18
    .1 .2100 .2500 0 .35 0
    .1 .0460 .0460 0 .1 0
    .1 .0460 .0460 0 -.1 0
    .1 .0460 .0230 -.08 -.605 0
    .1 .0230 .0230 0 -.606 0
    .1 .0230 .0460 .06 -.605 0];
theta_num=length(theta);
P=zeros(N_d,theta_num); %存放投影数据
rho=shep(:,1); %椭圆对应密度
ae=0.5*N*shep(:,2);%椭圆短半轴
be=0.5*N*shep(:,3);%椭圆长半轴
xe=0.5*N*shep(:,4);%椭圆中心x坐标
ye=0.5*N*shep(:,5);%椭圆中心y坐标
alpha=shep(:,6);%椭圆旋转角度
alpha=alpha*pi/180;
theta=theta*pi/180;%角度换算弧度
TT=-(N_d-1)/2:(N_d-1)/2;%探测器坐标
for k1=1:theta_num
    P_theta=zeros(1,N_d);
    for k2=1:max(size(xe))
        a =(ae(k2)*cos(theta(k1)-alpha(k2)))^2+(be(k2)*sin(theta(k1)-alpha(k2)))^2;
        temp=a-(TT-xe(k2)*cos(theta(k1))-ye(k2)*sin(theta(k1))).^2;
        ind=temp>0;%根号内值需为非负
        P_theta(ind)=P_theta(ind)+rho(k2)*(2*ae(k2)*be(k2)*sqrt(temp(ind)))./a;
    end
    P(:,k1)=P_theta;
end