如果您想在傅立叶域中进行不带滤波的逆 Radon 变换，则您列出的公式是一个中间结果.另一种方法是在空间域中使用卷积来完成整个滤波反投影算法，这在 40 年前可能会更快；你最终会重新推导你发布的公式.但是，我现在不推荐它，尤其是对于您的第一次重建；你应该先真正了解希尔伯特变换.

The formula you listed is an intermediate result if you wanted to do an inverse Radon transform without filtering in the Fourier domain. An alternative is to do the entire filtered back projection algorithm using convolution in the spatial domain, which might have been faster 40 years ago; you would eventually rederive the formula you posted. However, I wouldn't recommended it now, especially not for your first reconstruction; you should really understand the Hilbert transform first.

无论如何，这里有一些 Matlab 代码，它执行强制性的 Shepp-Logan 幻像滤波反投影重建.我将展示如何使用 Ram-Lak 过滤器进行自己的过滤.如果我真的有动力，我会用一些 interp2 命令和求和来替换 radon/iradon.

Anyway, here's some Matlab code which does the obligatory Shepp-Logan phantom filtered back projection reconstruction. I show how you can do your own filtering with the Ram-Lak filter. If I was really motivated, I would replace radon/iradon with some interp2 commands and summations.

phantomData=phantom();

phantomData=phantom();

N=size(phantomData,1);

theta = 0:179;
N_theta = length(theta);
[R,xp] = radon(phantomData,theta);

% make a Ram-Lak filter. it's just abs(f).
N1 = length(xp);
freqs=linspace(-1, 1, N1).';
myFilter = abs( freqs );
myFilter = repmat(myFilter, [1 N_theta]);

% do my own FT domain filtering
ft_R = fftshift(fft(R,[],1),1);
filteredProj = ft_R .* myFilter;
filteredProj = ifftshift(filteredProj,1);
ift_R = real(ifft(filteredProj,[],1));

% tell matlab to do inverse FBP without a filter
I1 = iradon(ift_R, theta, 'linear', 'none', 1.0, N);

subplot(1,3,1);imagesc( real(I1) ); title('Manual filtering')
colormap(gray(256)); axis image; axis off

% for comparison, ask matlab to use their Ram-Lak filter implementation
I2 = iradon(R, theta, 'linear', 'Ram-Lak', 1.0, N);

subplot(1,3,2);imagesc( real(I2) ); title('Matlab filtering')
colormap(gray(256)); axis image; axis off

% for fun, redo the filtering wrong on purpose
% exclude high frequencies to create a low-resolution reconstruction
myFilter( myFilter > 0.1 ) = 0;
ift_R = real(ifft(ifftshift(ft_R .* myFilter,1),[],1));
I3 = iradon(ift_R, theta, 'linear', 'none', 1.0, N);

subplot(1,3,3);imagesc( real(I3) ); title('Low resolution filtering')
colormap(gray(256)); axis image; axis off