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

blue的Whale鲸鱼

已于 2024-06-07 17:08:03 修改

阅读量628

点赞数 11

分类专栏： CT图像重建算法文章标签：算法 matlab 重构

于 2024-06-07 10:32:02 首次发布

本文链接：https://blog.csdn.net/2201_76025250/article/details/139519904

版权

CT图像重建算法专栏收录该内容

4 篇文章 0 订阅

订阅专栏

联合代数重建算法(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)=backproj(ind2)./sumCol(ind2);
        F=F+lamda*delta;
        F(F<0)=0;
    end
    c=c+1;
end
F=reshape(F,N,N)';
figure(1);
imshow(I);xlabel('(a)180*180头模型图像');
figure(2);
imshow(F);xlabel('(b)SART算法重建图像');


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

%计算投影矩阵
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

blue的Whale鲸鱼

关注

11
点赞
踩
5

收藏

觉得还不错? 一键收藏
0
评论
CT图像重建算法——联合代数重建算法(SART)

不同的是，在一次迭代更新中SART算法使用的是某一投影角度下通过某一像素的所有射线.而SIRT算法用到的是通过该像素的所有射线。SART算法对ART算法做了一些改进：利用同一投影角度下通过像素的所有射线的误差来确定对该个像素进行校正值，而不是只考虑一条射线。SART的迭代公式为。联合代数重建算法(simultaneous ART,SART)是A.H.Andersen和A.C.Kak于1984年提出一种改进的迭代重建算法，只需要很少的选代次数就可以得到很好的重建质量和精度。
复制链接

扫一扫