我最近一直在玩断层重建算法。我已经有了很好的 FBP、ART、类似 SIRT/SART 的迭代方案的工作实现，甚至使用直线代数（慢！）。 这个问题与任何这些技术无关；“为什么有人会这样做，这里有一些 FBP 代码”形式的答案不是我想要的。
我想用这个程序做的下一件事是“完成集合”并实现所谓的“傅里叶重建方法”。我对此的理解基本上是您将一维 FFT 应用于正弦图“曝光”，将它们排列为二维傅里叶空间中的径向“轮辐”（这是直接从中心切片定理得出的有用的事情） ，从这些点插入到该二维空间中的规则网格，然后应该可以进行傅里叶逆变换以恢复原始扫描目标。
听起来很简单，但我没有太多运气得到任何看起来像原始目标的重建。
下面的 Python (numpy/SciPy/Matplotlib) 代码是关于我想要做的最简洁的表达。运行时，它显示以下内容：
图1：目标
图2：目标的正弦图
图 3：FFT-ed 正弦图行
图 4：顶行是从傅里叶域正弦图行内插的 2D FFT 空间；底行是（用于比较目的）目标的直接 2D FFT。这是我开始怀疑的地方。从正弦图 FFT 插值的图看起来类似于直接对目标进行 2D-FFT 处理的图......但又有所不同。
图 5：图 4 的傅里叶逆变换。我希望这比实际更容易识别为目标。
任何想法我做错了什么？不确定我对傅立叶方法重建的理解是否存在根本缺陷，或者我的代码中只是存在一些错误。
import math
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate
import scipy.fftpack
import scipy.ndimage.interpolation
S=256 # Size of target, and resolution of Fourier space
A=359 # Number of sinogram exposures
# Construct a simple test target
target=np.zeros((S,S))
target[S/3:2*S/3,S/3:2*S/3]=0.5
target[120:136,100:116]=1.0
plt.figure()
plt.title("Target")
plt.imshow(target)
# Project the sinogram
sinogram=np.array([
np.sum(
scipy.ndimage.interpolation.rotate(
target,a,order=1,reshape=False,mode='constant',cval=0.0
)
,axis=1
) for a in xrange(A)
])
plt.figure()
plt.title("Sinogram")
plt.imshow(sinogram)
# Fourier transform the rows of the sinogram
sinogram_fft_rows=scipy.fftpack.fftshift(
scipy.fftpack.fft(sinogram),
axes=1
)
plt.figure()
plt.subplot(121)
plt.title("Sinogram rows FFT (real)")
plt.imshow(np.real(np.real(sinogram_fft_rows)),vmin=-50,vmax=50)
plt.subplot(122)
plt.title("Sinogram rows FFT (imag)")
plt.imshow(np.real(np.imag(sinogram_fft_rows)),vmin=-50,vmax=50)
# Coordinates of sinogram FFT-ed rows' samples in 2D FFT space
a=(2.0*math.pi/A)*np.arange(A)
r=np.arange(S)-S/2
r,a=np.meshgrid(r,a)
r=r.flatten()
a=a.flatten()
srcx=(S/2)+r*np.cos(a)
srcy=(S/2)+r*np.sin(a)
# Coordinates of regular grid in 2D FFT space
dstx,dsty=np.meshgrid(np.arange(S),np.arange(S))
dstx=dstx.flatten()
dsty=dsty.flatten()
# Let the central slice theorem work its magic!
# Interpolate the 2D Fourier space grid from the transformed sinogram rows
fft2_real=scipy.interpolate.griddata(
(srcy,srcx),
np.real(sinogram_fft_rows).flatten(),
(dsty,dstx),
method='cubic',
fill_value=0.0
).reshape((S,S))
fft2_imag=scipy.interpolate.griddata(
(srcy,srcx),
np.imag(sinogram_fft_rows).flatten(),
(dsty,dstx),
method='cubic',
fill_value=0.0
).reshape((S,S))
plt.figure()
plt.suptitle("FFT2 space")
plt.subplot(221)
plt.title("Recon (real)")
plt.imshow(fft2_real,vmin=-10,vmax=10)
plt.subplot(222)
plt.title("Recon (imag)")
plt.imshow(fft2_imag,vmin=-10,vmax=10)
# Show 2D FFT of target, just for comparison
expected_fft2=scipy.fftpack.fftshift(scipy.fftpack.fft2(target))
plt.subplot(223)
plt.title("Expected (real)")
plt.imshow(np.real(expected_fft2),vmin=-10,vmax=10)
plt.subplot(224)
plt.title("Expected (imag)")
plt.imshow(np.imag(expected_fft2),vmin=-10,vmax=10)
# Transform from 2D Fourier space back to a reconstruction of the target
fft2=scipy.fftpack.ifftshift(fft2_real+1.0j*fft2_imag)
recon=np.real(scipy.fftpack.ifft2(fft2))
plt.figure()
plt.title("Reconstruction")
plt.imshow(recon,vmin=0.0,vmax=1.0)
plt.show()
