matlab_医学CT重建 ART，SART算法

一，代码下载

matlab_医学CT重建ART，SART算法-自然语言处理文档类资源-CSDN下载

二、ART算法

1、基本思想

ART迭代算法的基本思想是先将连续的图像离散化，再采用CT成像的离散模型重建图像。其给定初始图像，先通过正投影得到投影图像，然后计算当前投影与实际测量投影之间的误差用以估计当前图像的修正值，这个修正值是对每一条射线逐条迭代修正并分配到射线穿过的像素上，再进行反投影和累加等处理。

2、算法实现步骤

• （1）给未知图像赋初值，例如：

fjk=0,k=0,j=1,2,⋯,N

• （2）第 i 条射线，计算其估计投影值

pi∗=∑n=1Nwinfn(k)

• （3）计算实际投影与估计投影的误差

Δi=pi−pi∗

• （4）对第 j 个像素点的修正值

Cj=Δi∑n=1Nwin2wij

• （5）对第 j 个像素点的值进行修正

fj(k+1)=fj(k)+λCj

• （6）将上一轮的结果作为初值，重复（2）~（5）的过程，直到达到收敛要求或指定的迭代次数。

因此对i条射线，ART的迭代公式为：

fj(k+1)=fj(k)+λpi−∑n=1Nwinfn(k)∑n=1Nwin2wij,j=1,2,⋯,N

其中， k 为迭代次数， i 为射线序号 (1<i<M) ， j 为像素序号 （）（1≤j≤N） , λ 为松弛因子 (0<λ<2) 。逐射线更新，每一轮的迭代都是以上一轮迭代的结果作为初始值，直到满足一定的收敛条件。

3、优缺点算法

从算法实现的步骤中可以看出，所谓的修正实际上只是针对射线所穿过的那些像素逬行的，而对于射线没有穿过的像素， wij=0 ，因而 Cj=0 ,相当于没有做出修正，因此其收敛速度慢，运算时间长。

ART 算法性能的影响因素包括了投影数据的读取方式和松弛因子的选择。由于每次迭代只用到一条射线投影，算法容易被测量数据中含有的噪声，该算法的这一缺点一定程度上限制了其应用和发展。

另一方面，对于ART算法，不同角度投影之间的先后次序对重建效果有一定影响。选择接近正交的方程，其收敛所需的迭代次数越小，收敛速度越快。

三、SART算法

1、基本思想

SART算法是A.H.Anderson和A.C.Kak于1984年提出一种改进的迭代重建算法，只需要很少的迭代次数就可以得到很好的重建质量和精度。SART算法对ART算法做了一些改进：利用同一投影角度下通过像素的所有射线的误差来确定对该个像素进行校正值，而不是只考虑一条射线。其效果相当于是对 ART 算法中的误差进行了平滑处理，从而降低重建结果对测量误差的敏感度。

2、算法实现步骤

• （1）对第 i 条射线，计算估计投影值

pi∗=∑n=1Nwinfn(k)

• （2）计算实际投影与估计投影的误差

Δi=pi−pi∗=pi−∑n=1Nwinfn(k)

• （3）反投影值

Cj=∑i=1Mpi−∑n=1Nwinfn(k)∑n=1Nwinwij

• （4）对第 j 个像素点的值进行修正

fj(k+1)=fj(k)+λCj∑i=1Mwij

• （5）将上一轮的结果作为初值，重复（1）~（4）的过程，直到达到收敛要求或指定的迭代次数。

因此，SART算法的迭代公式为：

fj(k+1)=fj(k)+λ∑i=1Mpi−∑n=1Nwinfn(k)∑n=1Nwinwij∑i=1Mwij

其中， λ 是松弛因子， k 是迭代次数。

3、优缺点分析

跟ART算法相比较，SART算法不是按照逐条射线对图像像素进行更新，而是在计算完一个特定投影角度的整个投影之后再进行更新。因此，SART算法的计算效率更高。且在每个投影角度下对像素点进行修正时，ART算法中只考虑了一条射线，与之相比，SART算法则是利用了通过像素点的所有射线，稳定性和并行性都更高。

SIRT算法和SART算法都能控制测量误差和干扰因素在迭代重建中的传播。虽然SIRT算法会使得像素点修正在一定程度上有所改善，但在计算量方面会比 SART算法大，因此，在实际应用中广泛使用的算法还是SART算法。

另一方面，对于SART算法，不同角度投影之间的先后次序对重建效果也有一定影响。选择接近正交的方程，其收敛所需的迭代次数越小，收敛速度越快。

四，部分源码

1，ART算法部分源码

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');%旋转采样
        ......

2，SART算法部分源码

clear all
close all
tic
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;
figure;
N=10;
angle=0:1:179;   
[R,xp]  = radon(I, angle);
R(200,1:180)=100;
[nx,ny] = size(I);
nb = size(xp,1);
na= size(angle,2);
W = W(nx,ny,nb,na);
r_sart_c = sart(W, ones(nx*ny,1), R(:),N,nx);
......

五，运行结果

1，ART结果

2，SART结果

备注：算法介绍部分摘自互联网，如有侵权，请联系删除