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

程序如下：

clear all

clc

I=imread('up4-Amp.png')

OutImg=I

R=I(:,:,1)

G=I(:,:,2)

B=I(:,:,3)

R=medfilt2(R,[3,3])

G=medfilt2(G,[3,3])

B=medfilt2(B,[3,3])

I1=cat(3,R,G,B)  % 对彩色图像R,G，B三个通道分别进行3×3模板的中值滤波 cat函数用于连接两个矩阵或数组

R=filter2(fspecial('average',3),R)/255

G=filter2(fspecial('average',3),G)/255

B=filter2(fspecial('average',3),B)/255

I2= cat(3,R,G,B)  %对彩色图像R,G，B三个通道分别进行3×3模板的均值滤波

figure,imshow(I)

title('原图')

figure,

imshow(I1)

title('中值滤波')

figure,imshow(I2)

title('均值滤波')

扩展资料：

注意事项

1、在频域滤波，由于是点乘，所以滤波模板矩阵和图像矩阵必须尺寸一样。

2、因为尺寸一样，它们的原点必须要对齐。

3、因在进行离散傅里叶变换后，在频域点乘，相当于在时域卷积，但是这个时候实际上是对时樱仿域周期矩阵进行卷积。直接脊粗纤在时域卷积，matlab默认是在边界补0。

4、Matlab freqz2（）这个函数可以自动得到一个指定尺寸的，对应于时域的频域模板。

5、图像经过傅里叶变换后，它的原点在左上角。而模板经过freqz2后，原点在凳磨中心，所以只要平移其中的一个就好了。

6、在对原图像进行傅里叶变换之前，按照一定规则补0就好了。

程序是下面这样，但只能处理长宽一样的方形图像，灰度和彩祥宽谨色图像都可，要用其他图像只需把Lena.bmp改为其他图像，但图像要保存在m文件所在路径下。下面巧耐还有一个运行后的图像（之一），希望能对你有所帮助。

clearclc

%%%%%%%%%%测试图像只能是方形图像，长宽像素一样。

f=imread('Lena.bmp')%%读取图像数据，图像只能保存在m文件所在的路径下

d=size(f)

if length(d)>2

f=rgb2gray((f))%%%%%%%%如果是彩色图像则转化为灰度图

end

T=d(1)

SUB_T=T/2

% 2.进行二维小波分解

l=wfilters('db10','l') % db10（消失矩为10)低通分解滤波器冲击响应（长度为20）

L=T-length(l)

l_zeros=[l,zeros(1,L)] % 矩阵行数与输入图像一致，为2的整数幂

h=wfilters('db10','h') % db10（消失矩为10)高通分解滤波器冲击响应（长度为20）

h_zeros=[h,zeros(1,L)] % 矩阵行数与输入图像一致，为2的整数幂

for i=1:T % 列变换

row(1:SUB_T,i)=dyaddown( ifft( fft(l_zeros).*fft(f(:,i)') ) ).' % 圆周卷积<->FFT

row(SUB_T+1:T,i)=dyaddown( ifft( fft(h_zeros).*fft(f(:,i)') ) ).' % 圆周卷积<->FFT

end

for j=1:T % 行变换

line(j,1:SUB_T)=dyaddown( ifft( fft(l_zeros).*fft(row(j,:)) ) ) % 圆周卷积<->FFT

line(j,SUB_T+1:T)=dyaddown( ifft( fft(h_zeros).*fft(row(j,:)) ) ) % 圆周谨基卷积<->FFT

end

decompose_pic=line % 分解矩阵

% 图像分为四块

lt_pic=decompose_pic(1:SUB_T,1:SUB_T) % 在矩阵左上方为低频分量--fi(x)*fi(y)

rt_pic=decompose_pic(1:SUB_T,SUB_T+1:T) % 矩阵右上为--fi(x)*psi(y)

lb_pic=decompose_pic(SUB_T+1:T,1:SUB_T) % 矩阵左下为--psi(x)*fi(y)

rb_pic=decompose_pic(SUB_T+1:T,SUB_T+1:T) % 右下方为高频分量--psi(x)*psi(y)

% 3.分解结果显示

figure(1)

subplot(2,1,1)

imshow(f,[]) % 原始图像

title('original pic')

subplot(2,1,2)

image(abs(decompose_pic)) % 分解后图像

title('decomposed pic')

figure(2)

% colormap(map)

subplot(2,2,1)

imshow(abs(lt_pic),[]) % 左上方为低频分量--fi(x)*fi(y)

title('\Phi(x)*\Phi(y)')

subplot(2,2,2)

imshow(abs(rt_pic),[]) % 矩阵右上为--fi(x)*psi(y)

title('\Phi(x)*\Psi(y)')

subplot(2,2,3)

imshow(abs(lb_pic),[]) % 矩阵左下为--psi(x)*fi(y)

title('\Psi(x)*\Phi(y)')

subplot(2,2,4)

imshow(abs(rb_pic),[]) % 右下方为高频分量--psi(x)*psi(y)

title('\Psi(x)*\Psi(y)')

% 5.重构源图像及结果显示

% construct_pic=decompose_matrix'*decompose_pic*decompose_matrix

l_re=l_zeros(end:-1:1) % 重构低通滤波

l_r=circshift(l_re',1)' % 位置调整

h_re=h_zeros(end:-1:1) % 重构高通滤波

h_r=circshift(h_re',1)' % 位置调整

top_pic=[lt_pic,rt_pic] % 图像上半部分

t=0

for i=1:T % 行插值低频

if (mod(i,2)==0)

topll(i,:)=top_pic(t,:)% 偶数行保持

else

t=t+1

topll(i,:)=zeros(1,T) % 奇数行为零

end

end

for i=1:T % 列变换

topcl_re(:,i)=ifft( fft(l_r).*fft(topll(:,i)') )' % 圆周卷积<->FFT

end

bottom_pic=[lb_pic,rb_pic] % 图像下半部分

t=0

for i=1:T % 行插值高频

if (mod(i,2)==0)

bottomlh(i,:)=bottom_pic(t,:) % 偶数行保持

else

bottomlh(i,:)=zeros(1,T) % 奇数行为零

t=t+1

end

end

这个只是一级分解，matlab自带的函数可以实现多级分解，级数由编程者自己确定。

是的，是一样的。

---