滤波反投影（Filtered Back-Projection, FBP）是一种广泛应用于计算机断层扫描（CT）等影像技术中的重建算法。该算法基于反投影和滤波两个关键步骤，从一系列的投影数据中重建出物体的内部结构。
投影数据采集：首先，通过不同的角度对物体进行射线投影，收集到一系列的投影数据。这些数据反映了物体在不同方向上的密度分布信息。
滤波处理：投影数据在反投影之前会经过一个滤波过程。这个滤波函数通常选择为斜坡滤波器（Ramp Filter），用于减少低频噪声并增强图像的边缘信息，从而提高重建图像的清晰度。
反投影：经过滤波处理后的投影数据被反投影回物体空间。这个过程是通过计算每个投影角度下射线穿过物体的路径，并将滤波后的投影数据按照这些路径均匀地分布回去，从而得到一个初步的重建图像。
图像重建：通过叠加所有角度的反投影数据，最终得到物体的三维重建图像。
滤波反投影定理: 设 f(x,y) 是一个具有很好性质的函数， 它的拉东变换是
,
那么
证明 函数f的傅里叶变换形式为，
其逆变换为
。
第二个等号我们使用了变量变换
它的雅可比行列式定义为，
因此，
。
下一步推导依赖于投影-切片定理（Projection-Slice Theorem），它是CT图像重建中的一个关键定理。首先， 根据投影-切片定理，我们知道函数 f(x, y) 在某个角度 θ上的投影 g(t, θ) 的一维傅里叶变换等于f(x, y)的二维傅里叶变换在相应方向上的切片(后面会仔细解释这句话)。
变量替换：
在逆傅里叶变换公式中，我们将
表示为极坐标形式，即
，
其中 v 是径向频率，θ 是方向角度。这样做是为了将二维频率空间转换为极坐标形式，便于与投影数据对应起来。
代入投影数据: 由于投影-切片定理，我们知道在某个角度 θ 上，（即二维傅里叶变换在角度 θ 上的切片）与投影数据 g(t, θ) 的一维傅里叶变换相对应。因此，我们可以将
替换为投影数据的一维傅里叶变换形式。
在投影-切片定理的中，
表示二维函数 f(x, y) 的二维傅里叶变换在极坐标下的表示。根据投影-切片定理，对于每个固定的角度 θ，函数 f(x, y) 沿着这个角度的投影（Radon 变换）的一维傅里叶变换，等于
。
具体来说，如果我们定义 g(t, θ) 为函数 f(x, y) 在角度 θ 上的投影，即 Radon 变换，那么 g(t, θ) 的一维傅里叶变换（关于变量 t）就等于
。
用数学符号表示就是：
其中
表示对 g(t,θ) 关于 t 进行一维傅里叶变换
当然, 这一点如果觉得不自然, 可以看看 g(t,θ)的定义,
,
g(t,θ) 关于 t 进行一维傅里叶变换是
然后再用一个坐标变换
并且注意到 这个坐标变换的Jacobian 是1, 于是, 我们上面的积分换成了
推导这个式子用到了 从 导出的
.
再注意到 , 上式的右端的确就是 f(x,y)的傅里叶变换
.
现在，回到我们的推导过程中，由于我们需要从投影数据中重建原始图像 f(x, y)，我们可以利用投影-切片定理，将二维傅里叶变换 在每个方向
上的切片，替换为对应投影数据 g(t,θ) 的一维傅里叶变换。
因此，在逆傅里叶变换的公式中，当我们遇到
时，我们可以根据投影-切片定理，将其替换为 g(t,θ) 的一维傅里叶变换。
这样做的好处是，我们可以通过测量得到的投影数据（即 Radon 变换的结果），来间接地获取二维傅里叶变换的信息，进而通过逆傅里叶变换重建出原始的图像 f(x, y)。
逆变换中的积分: 在逆傅里叶变换中，我们需要对所有的频率和方向进行积分，以重建原始函数 f(x, y)。
因此，我们有两个积分：一个是对频率 v 的积分，另一个是对方向 θ 的积分。
上面的被积函数中的 v 来自于从笛卡尔坐标(v_x, v_y)转换到极坐标 (v, θ)时的雅可比行列式的值是 v（即径向频率）。 这里的 将要换成
。
通过上述推导，我们得到了公式：
注意，在实际应用中，v 通常会被替换为相应的滤波函数，以改善图像重建的质量。此外，公式中的
项表示对投影数据进行一维傅里叶变换（这是在频率域中进行的操作），而
项则表示逆傅里叶变换，将处理后的数据从频率域转换回空间域。
上面关于v积分范围从0到∞, θ积分范围从0到2π, 这和v积分范围从-∞到∞, θ积分范围从0到π 是等价的，这里给出一些解释。
首先，对于v的积分范围，从−∞到∞是二维傅里叶变换的标准形式。然而，在滤波后投影定理中，我们通常只考虑正频率（因为傅里叶变换的对称性），因此积分范围可以限制为从0到∞。这是因为在极坐标下，负频率可以通过旋转角度θ来覆盖，即对于每个负频率(−v,θ)，都存在一个正频率(v,θ+π)与之对应。
总结来说，两种积分范围是等价的，但选择哪种范围取决于具体的实现和上下文。在滤波后投影定理中，通常选择v从0到∞和θ从0到2π，因为这样更直接地与投影数据和滤波操作相对应但也v从−∞到∞和θ从0到π的范围确实是一种常见的表示方法。
我参考 GitHub - YuZitong/Radon-Transform-and-Inverse-Radon-Transform: Radon Transform and Inverse Radon Transform
把这儿一个代码移植到北太天元上
% 清除工作空间中的所有变量和关闭所有图形窗口
clear;
close all;
% 创建一个简单的矩形图像
imageSize = 12; % 图像的大小
simpleImage = zeros(imageSize, imageSize); % 创建一个全零矩阵
% 在图像中心绘制一个白色矩形区域
rectStartX = round(imageSize / 4);
rectEndX = round(3 * imageSize / 4);
rectStartY = round(imageSize / 4);
rectEndY = round(3 * imageSize / 4);
simpleImage(rectStartX:rectEndX, rectStartY:rectEndY) = 1;
% 图像大小
[r,c] = size(simpleImage);
% 旋转中心，即图像的中心点
xOrigin = floor(r/2);
yOrigin = floor(c/2);
% 投影角度
num_projection = 180; % 投影数量
delta_theta = 180/num_projection; % 角度增量
theta = linspace(0,180-delta_theta,num_projection); % 生成投影角度数组（不包括180度，以避免重复）
% 投影的大小
pFirst = floor(sqrt((r-xOrigin)^2+(c-yOrigin)^2))+1; % 计算最大投影长度并加1（为了包含边界）
rSize = 2*pFirst+1; % 投影数组的大小（考虑到正负偏移）
% 初始化投影矩阵
R = zeros(rSize, length(theta));
% 循环每个投影角度
for ag = theta
numAngle = ag/180*pi; % 将角度转换为弧度
SinA = sin(numAngle);
CosA = cos(numAngle);
% 构建cos和sin表以节省计算时间（但这里的方法可能不是最高效的）
xcos = zeros(2*c,1);
ysin = zeros(2*r,1);
% 计算每行对应的sin值
for i = 1:r
y = yOrigin - i;
ysin(2*i-1) = (y-0.25)*SinA; % 上半像素
ysin(2*i) = (y+0.25)*SinA; % 下半像素
end
% 计算每列对应的cos值
for j = 1:c
x = j - xOrigin;
xcos(2*j-1) = (x-0.25)*CosA; % 左半像素
xcos(2*j) = (x+0.25)*CosA; % 右半像素
end
% 遍历图像的每个像素
for j = 1:c % 列（x方向）
for i = 1:r % 行（y方向）
if simpleImage(i,j) ~= 0
p = double(simpleImage(i,j))*0.25; % 将一个像素拆分为4个子像素以提高精度
% 计算子像素在投影中的位置，并更新投影矩阵
tho = xcos(2*j-1)+ysin(2*i-1)+pFirst; R = contribute(R,tho,ag,delta_theta,p); % 左上
tho = xcos(2*j)+ysin(2*i-1)+pFirst; R = contribute(R,tho,ag,delta_theta,p); % 右上
tho = xcos(2*j-1)+ysin(2*i)+pFirst; R = contribute(R,tho,ag,delta_theta,p); % 左下
tho = xcos(2*j)+ysin(2*i)+pFirst; R = contribute(R,tho,ag,delta_theta,p); % 右下
end
end
end
end
% 显示投影（即sinogram）
imagesc(R);title('Radon变换结果'); colorbar;
% 定义辅助函数：计算单个子像素对投影的贡献
function R = contribute(R, tho, ag, delta_theta, p)
% 如果tho大于0，则取整数值，否则取整数值减1（为了处理负索引）
if tho > 0
Rr = floor(tho);
else
Rr = floor(tho)-1;
end
delta = tho - Rr; % 计算与整数值的偏差
% 按比例将贡献分配到最近的两个投影箱中
R(Rr+1, round(ag/delta_theta)+1) = R(Rr+1, round(ag/delta_theta)+1) + p*(1-delta); % 注意：这里使用round来确保索引为整数
R(Rr+2, round(ag/delta_theta)+1) = R(Rr+2, round(ag/delta_theta)+1) + p*delta;
end
可以把上面的radon 变换的脚本改造成子函数 my_radon, 然后来验证这篇文章中证明的一些结论.
% 清除工作空间中的所有变量和关闭所有图形窗口
clear;
close all;
% 1. 定义函数 f(x,y)，这里以高斯函数为例
[x, y] = meshgrid(-1.5:0.1:1.5, -1.5:0.1:1.5); % 定义网格
f = double( x.^2+y.^2 <= 1); % 单位圆
% 2. 计算 f(x,y) 的二维 Fourier 变换
F = fft2(f);
Fshift = fftshift(F); % 将零频分量移到频谱中心
% 3. 计算 f(x,y) 的 Radon 变换（拉东变换）
theta = 0:180; % 定义投影角度范围
R = my_radon(f); % 使用my_radon函数计算Radon变换
% 4. 验证投影切片定理
% 选择一个角度（例如 theta = 0）
theta_slice = 0; % 选择切片角度
idx_slice = find(theta == theta_slice, 1, 'first'); % 找到对应索引
% 计算二维Fourier变换在选定角度的切片
F_slice = Fshift(round(size(Fshift, 1) / 2), :); % 假设零频分量在中心行
% 计算Radon变换的一维Fourier变换
R_theta = R(:, idx_slice); % 获取选定角度的Radon变换数据
R_theta_fft = fft(R_theta); % 计算一维Fourier变换
R_theta_fftshift = fftshift(R_theta_fft); % 将零频分量移到频谱中心
% 可视化比较
figure;
subplot(2, 2, 1); imagesc(abs(Fshift)); title('2D Fourier Transform of f(x,y)'); colormap gray;
axis square; colorbar;
subplot(2, 2, 2); imagesc(abs(R)); title('Radon Transform of f(x,y)'); colormap gray;
axis square; colorbar;
subplot(2, 2, 3); plot(abs(F_slice)); title('Slice of 2D Fourier Transform'); grid on;
xlabel('Frequency'); ylabel('Magnitude');
subplot(2, 2, 4); plot(abs(R_theta_fftshift)); title('1D Fourier Transform of Radon Slice'); grid on;
xlabel('Frequency'); ylabel('Magnitude');
function R = my_radon(simpleImage)
% 图像大小
[r,c] = size(simpleImage);
% 旋转中心，即图像的中心点
xOrigin = floor(r/2);
yOrigin = floor(c/2);
% 投影角度
num_projection = 180; % 投影数量
delta_theta = 180/num_projection; % 角度增量
theta = linspace(0,180-delta_theta,num_projection); % 生成投影角度数组（不包括180度，以避免重复）
% 投影的大小
pFirst = floor(sqrt((r-xOrigin)^2+(c-yOrigin)^2))+1; % 计算最大投影长度并加1（为了包含边界）
rSize = 2*pFirst+1; % 投影数组的大小（考虑到正负偏移）
% 初始化投影矩阵
R = zeros(rSize, length(theta));
% 循环每个投影角度
for ag = theta
numAngle = ag/180*pi; % 将角度转换为弧度
SinA = sin(numAngle);
CosA = cos(numAngle);
% 构建cos和sin表以节省计算时间（但这里的方法可能不是最高效的）
xcos = zeros(2*c,1);
ysin = zeros(2*r,1);
% 计算每行对应的sin值
for i = 1:r
y = yOrigin - i;
ysin(2*i-1) = (y-0.25)*SinA; % 上半像素
ysin(2*i) = (y+0.25)*SinA; % 下半像素
end
% 计算每列对应的cos值
for j = 1:c
x = j - xOrigin;
xcos(2*j-1) = (x-0.25)*CosA; % 左半像素
xcos(2*j) = (x+0.25)*CosA; % 右半像素
end
% 遍历图像的每个像素
for j = 1:c % 列（x方向）
for i = 1:r % 行（y方向）
if simpleImage(i,j) ~= 0
p = double(simpleImage(i,j))*0.25; % 将一个像素拆分为4个子像素以提高精度
% 计算子像素在投影中的位置，并更新投影矩阵
tho = xcos(2*j-1)+ysin(2*i-1)+pFirst; R = contribute(R,tho,ag,delta_theta,p); % 左上
tho = xcos(2*j)+ysin(2*i-1)+pFirst; R = contribute(R,tho,ag,delta_theta,p); % 右上
tho = xcos(2*j-1)+ysin(2*i)+pFirst; R = contribute(R,tho,ag,delta_theta,p); % 左下
tho = xcos(2*j)+ysin(2*i)+pFirst; R = contribute(R,tho,ag,delta_theta,p); % 右下
end
end
end
end
end
% 定义辅助函数：计算单个子像素对投影的贡献
function R = contribute(R, tho, ag, delta_theta, p)
% 如果tho大于0，则取整数值，否则取整数值减1（为了处理负索引）
if tho > 0
Rr = floor(tho);
else
Rr = floor(tho)-1;
end
delta = tho - Rr; % 计算与整数值的偏差
% 按比例将贡献分配到最近的两个投影箱中
R(Rr+1, round(ag/delta_theta)+1) = R(Rr+1, round(ag/delta_theta)+1) + p*(1-delta); % 注意：这里使用round来确保索引为整数
R(Rr+2, round(ag/delta_theta)+1) = R(Rr+2, round(ag/delta_theta)+1) + p*delta;
end
