求图像重建反投影、滤波反投影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('滤波反投迅扮影法')

把采集到的图象用仿射变换配准，

为了加快运行速度可以先敏山猜进行展开。

配准这一步可以在空间域，

也可在频率域进行

然后按配准结果将这些图象插合成一幅图象，

再用最小二乘法求解线性方程组即可。

注意，

最好使用超松弛迭代法求解，

但是遇到0的时候结果可能有较大唯亏出入，

解决办法中的一种是图象矩阵所有元素全部加上1，

计算完成后再全部减去1，

然后再512级灰度量化

这是最简单的重构方法，

没有考虑图象的模糊效应。

此外，如果有矩阵维度问题，

有两种解决办法，

一是将插合图象变成正方形图象，

一是将各插合行，列按权值累加，

反向桥型映射，

后一种速度快些，

也不必直接求解方程，

但是不具有通用性。

---