求图像重建反投影、滤波反投影matlab仿真程序

%用phantom函数可以获得仿体图像；

%用randon可获得不同角度的一维投影；

clear all

P = phantom('Modified Shepp-Logan',256)

R=radon(P)

figureimshow(R,[])

figure

imshow(P,[])title('仿体图')

%直接反投影法

l = pow2(nextpow2(size(R,1))-1)%重构图像的大小

P_1 = zeros(l,l)%用于存放重构后的图像

for i=1 : size(R,2)

    tmp = imrotate( repmat(R(:,i),1,size(R,1)),i-1,'bilinear' )

    tmp = tmp(floor(size(tmp,1)/2-l/2)+1:floor(size(tmp,1)/2+l/2),floor(size(tmp,2)/2-l/2)+1:floor(size(tmp,2)/2+l/2))

    P_1=P_1+tmp

end

P_1=P_1/size(R,2)

P_1=rot90(P_1)

figureimshow(P_1,[])title('直接反投影法')

%滤波反投影法

N=180

%滤波

H=size(R,1)

h=zeros((H*2-1),1)

for i=0:H-1

    if i==0

        h(H-i)=1/4

    elseif rem(i,2)==0

        h(H-i)=0

  姿御      h(H+i)=0

    else

        h(H-i)=-1/(i*pi)^2

        h(H+i)=-1/(i*pi)^2

  迹昌岩  end

end

x=zeros(H,N)

for i=1:N

    s=R(:,i)

    xx=conv(s',h')

    x(:,i)=xx(H:2*H-1)

end

%反投影

P_3=zeros(l,l)

for i=1:l

    for j=1:l

        for k=1:180

            theta=k/180*pi

            t=(j-l/2-0.5)*cos(theta)+(l/2+0.5-i)*sin(theta)+(H+1)/2

            t1=floor(t)

            t2=floor(t+1) 

            P_3(i,j)=P_3(i,j)+(t2-t)*x(t1,k)+(t-t1)*x(t2,k)

        end

    end

end

P_3=pi/N*P_3

figureimshow(P_3,[])title('滤波反投迅扮影法')

iradon函数是基于R-L滤波器的滤波反投影法冲启。实现重建图像的过程如下：

1、把投影矩阵R转饥判型换到频域，生成fft（R）；

2、fft（R）和滤波函数H相乘，得到滤波后的频域投影矩阵fft（R）*H；

3、把fft（R）*H 转换到空域，得到空域中的滤波后的投影矩阵R'=ifft（fft（R）*H ）

4、R'插值后，得到处理好的投影矩阵R''

5、直接反投影得到重建图像矩阵I。

滤波函数H如果选择“none”，则没有滤波，选择“ram-lak”，则把R-L滤波函数的傅里叶函数，和频域中每烂猜个角度的投影相乘，实现滤波。选择其他，则R-L滤波函数的傅里叶函数与设定的函数相乘，再和频域中每个角度的投影相乘。

滤波反投影法与迭代重建算隐液法的优缺点比较：

确定迭代变量。在可以用迭代算法解决的问题中，至少存在一个直接拆此或间接地不断由旧值递推出新值的变量，这个变量就是迭代变量。

建立迭代关系式。所谓迭代关系式，指如何从变量的前一个值推出其下一个值的公式(或关系)。迭代关系式的建立是解决迭代问题的关键，通常可以使用递推或倒推的方法来完成。

过程控制

在什么时候结束迭代过程?这是编写迭代程序必须考虑的问题。不能让迭代过程无休止地重复执行下去。迭代过程的控制通常可分为两种情况：一种是所需的迭代次数是个确定的值，可以计算出来；另一种是所需的迭代次数旅携迅无法确定。对于前一种情况，可以构建一个固定次数的循环来实现对迭代过程的控制；对于后一种情况，需要进一步分析出用来结束迭代过程的条件。

---