关于投影与反投影的matlab代码

关于投影与反投影的在matlab代码中的运用，经供参考，谢谢。

clear all;
close all;

%% 获得原始图像、投影数据、以及直接反投影后的断层图像
p_N = 256; 
theta_N = 180; 
pad_N = 1024;

P = phantom（p_N）;
theta = 0:（theta_N-1）;
[Rxp] = radon（Ptheta）;
[I H] = iradon（Rtheta‘linear‘‘none‘）;                % iradon变换，不滤波，则I就是直接反投影后的断层图像

row_pad = pad_N; 
row_I = size（I1）; 
row_start = （row_pad-row_I）/2;
fb = zeros（row_padrow_pad）;                            % 将I扩充成1024*1024
for i = 1:1:row_I
    for j = 1:1:row_I
        fb（i+row_startj+row_start） = I（ij）;           % 里面256*256是直接反投影后的数据，外面的都是0
    end
end

%% 获取二维锥形滤波器Frow， Frow = sqrt（u^2+v^2）
N = pad_N; 
uP = zeros（NN）; 
vP = zeros（NN）; 
Frow = zeros（NN）;
for i = 1:1:N
    for j = 1:1:N
        uP（ij） = （（-（N）/2）+（j）*1）*2*pi/N;
        vP（ij） = （（-（N）/2）+（i）*1）*2*pi/N;
        Frow（ij） = sqrt（uP（ij）^2+vP（ij）^2）;
    end
end

%% 求取断层图像fb的二维傅里叶变换Fb，根据Frow二维滤波后得出Fvu，Fvu = Fb*Frow
Fb = fft2（fb）;                                          % 二维傅里叶变换
Fb = ifftshift（Fb）;                                     % 直流分量移动到中央，中心是低频，向外是高频
Fvu = zeros（row_padrow_pad）;
for i=1:1:row_pad
		for j=1:1:row_pad
			if（i==（row_pad/2+1） || j==（row_pad/2+1） ）	
				Fvu（ij） = Fb（ij）;                     % 中央的直流分量，仍然保持未滤波前的值
			else
				Fvu（ij） = Fb（ij）*Frow（ij）;			% 二维滤波
			end
		end
end

%% 对滤波后的数据Fvu，求取二维傅里叶逆变换，得到fxy
fxy = （ifft2（（Fvu）））;                                   % 求取二维傅里叶逆变换
target = fxy;
crop = row_start; 
target = target（crop+1:end-cropcrop+1:end-crop）;
I_a = abs（target）;                                      % 复数的模值 
I_a = （I_a-min（I_a（:）））./（max（I_a（:））-min（I_a（:）））;     % 归一化为 0~1
Lg_I_a = log（1+I_a）;                                    % 为了好的显示效果，取对数

figure 
subplot（131） imshow（P[]） title（‘原始图像‘）;
subplot（132） imshow（I[]） title（‘直接反投影后的断层图像‘）;
subplot（133） imshow（Lg_I_a[]） title（‘先反投影后滤波恢复的图像‘）;

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       2107  2017-06-08 09:37  二维图像的投影和图像重建分析之傅里叶变换法\BPF.m

     文件       1729  2017-06-08 09:41  二维图像的投影和图像重建分析之傅里叶变换法\BPF_conv2.m

     文件       7090  2017-06-04 19:16  二维图像的投影和图像重建分析之傅里叶变换法\CBP\Backproject.c

     文件       8704  2017-06-04 20:35  二维图像的投影和图像重建分析之傅里叶变换法\CBP\Backproject.mexw64

     文件       2737  2017-06-08 09:44  二维图像的投影和图像重建分析之傅里叶变换法\CBP\CBP.m

     文件       1686  2017-06-05 10:56  二维图像的投影和图像重建分析之傅里叶变换法\CBP\DirectBackProject.m

     文件       3394  2017-06-08 09:38  二维图像的投影和图像重建分析之傅里叶变换法\DFR.m

     文件        766  2017-06-04 16:46  二维图像的投影和图像重建分析之傅里叶变换法\FBP.m

     文件        383  2017-06-08 19:45  二维图像的投影和图像重建分析之傅里叶变换法\说明.txt

     目录          0  2018-08-28 16:21  二维图像的投影和图像重建分析之傅里叶变换法\CBP

     目录          0  2018-08-28 16:21  二维图像的投影和图像重建分析之傅里叶变换法

----------- ---------  ---------- -----  ----

                28596                    11