用傅立叶方法进行层析成像重建的代码有什么问题？

最近，我一直在研究层析重建算法。我已经有了FBP，ART，类似于SIRT / SART的迭代方案甚至使用直线线性代数的慢速工作实现（慢！）。这个问题与这些技术无关。 “为什么有人会这样做，下面是一些FBP代码”这样的形式的答案并不是我想要的。

我想用该程序做的下一件事是“完成集”，并实施所谓的“傅立叶重建法”。我对此的理解基本上是对正弦图“曝光”应用1D FFT，将其作为径向“轮辐”布置在2D傅立叶空间中（这是非常有用的事情，直接根据中心切片定理可以得出），从这些点插入到2D空间中的规则网格中，然后应该可以进行傅立叶逆变换来恢复原始扫描目标。

听起来很简单，但是我没有太多运气好的话，看起来像是原始目标的重构。

下面的Python（numpy / SciPy / Matplotlib）代码是我想出的最简洁的表达式。运行时，它将显示以下内容：图1：目标

图2：目标的正弦图

图3：FFT编码的正弦图行

图4：最上面的行是从傅立叶域正弦图行内插的2D FFT空间；最下面一行是（出于比较目的）目标的直接2D FFT。这是我开始变得可疑的地方。从正弦图FFT插值的图看起来与通过直接对目标进行2D FFT绘制的图相似，但又有所不同。
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()

这是否等同于使用FFT计算Radon逆变换？

...因为这里有代码，应该放在中心的东西在边缘，应该放在边缘的东西在中心，例如不应有90度的相移？

链接的代码用于过滤后向投影（FBP）方法。它基于相同的中心切片数学，但从未明确尝试构建2D傅里叶域图像。您可以将FBP滤波器对低频的抑制看作是对中央部分“辐条”的较高密度的补偿。在我尝试实现的傅里叶重构方法中，这只是表现为要插值的点的密度更高。我自由地承认我正在尝试实施一种很少使用的技术，并且在文献中对该技术的覆盖范围很有限，

糟糕，是的，您是对的。这是C语言中的一个版本。我仔细浏览了一下并发布了一些内容。我待会儿再看。

#1 楼

好吧，我终于破解了。

诀窍基本上归结为将一些fftshift / ifftshift放置在正确的位置，因此2D傅里叶空间表示法并不会产生剧烈的振荡，并且注定无法精确内插。至少这是我认为已解决的问题。我对傅立叶理论的有限理解大多数都是基于连续积分公式，而且我总是会发现离散域和FFT有点古怪。，我必须归功于这种实现方式，至少要给我信心，这种重建算法可以在这种环境中合理紧凑地表达。 >
图1：一个新的更复杂的目标。

图2：目标的正弦图（确定，这是Radon变换）。

图3：正弦图的FFT行（在中心绘制DC）。

图4：FFT变换的正弦图到2D FFT空间（DC居中）。颜色是绝对值的函数。

图4a：放大2D FFT空间的中心只是为了更好地显示正弦图数据的径向性质。

图5：顶行：从径向排列的FFT-ed正弦图行中插入的2D FFT空间。最下面一行：只需对目标进行2D FFT即可获得的预期外观。

图5a：放大图5中子图的中心区域以显示这些外观在质量上非常吻合。

图6：酸性测试：内插FFT空间的逆2D FFT恢复了目标。尽管我们已经完成了所有准备工作，但Lena的外观仍然很好（可能是因为有足够的正弦图“辐条”可以相当密集地覆盖2D FFT平面；如果减少曝光角度数，事情会变得很有趣，所以这不再成立了））。

代码如下；在i7上的Debian / Wheezy的64位SciPy上以不到15s的时间绘制图。

import math
import matplotlib
import matplotlib.pyplot as plt
import numpy as np

import scipy.interpolate
import scipy.fftpack
import scipy.misc
import scipy.ndimage.interpolation

S=256 # Size of target, and resolution of Fourier space
N=259 # Number of sinogram exposures (odd number avoids redundant direct opposites)

V=100 # Range on fft plots

# Convenience function
def sqr(x): return x*x

# Return the angle of the i-th (of 0-to-N-1) sinogram exposure in radians.
def angle(i): return (math.pi*i)/N

# Prepare a target image
x,y=np.meshgrid(np.arange(S)-S/2,np.arange(S)-S/2)
mask=(sqr(x)+sqr(y)<=sqr(S/2-10))
target=np.where(
    mask,
    scipy.misc.imresize(
        scipy.misc.lena(),
        (S,S),
        interp='cubic'
        ),
    np.zeros((S,S))
    )/255.0

plt.figure()
plt.title("Target")
plt.imshow(target)
plt.gray()

# Project the sinogram (ie calculate Radon transform)
sinogram=np.array([
        np.sum(
            scipy.ndimage.interpolation.rotate(
                target,
                np.rad2deg(angle(i)), # NB rotate takes degrees argument
                order=3,
                reshape=False,
                mode='constant',
                cval=0.0
                )
            ,axis=0
            ) for i in xrange(N)
        ])

plt.figure()
plt.title("Sinogram")
plt.imshow(sinogram)
plt.jet()

# Fourier transform the rows of the sinogram, move the DC component to the row's centre
sinogram_fft_rows=scipy.fftpack.fftshift(
    scipy.fftpack.fft(
        scipy.fftpack.ifftshift(
            sinogram,
            axes=1
            )
        ),
    axes=1
    )

plt.figure()
plt.subplot(121)
plt.title("Sinogram rows FFT (real)")
plt.imshow(np.real(sinogram_fft_rows),vmin=-V,vmax=V)
plt.subplot(122)
plt.title("Sinogram rows FFT (imag)")
plt.imshow(np.imag(sinogram_fft_rows),vmin=-V,vmax=V)

# Coordinates of sinogram FFT-ed rows' samples in 2D FFT space
a=np.array([angle(i) for i in xrange(N)])
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()

plt.figure()
plt.title("Sinogram samples in 2D FFT (abs)")
plt.scatter(
    srcx,
    srcy,
    c=np.absolute(sinogram_fft_rows.flatten()),
    marker='.',
    edgecolor='none',
    vmin=-V,
    vmax=V
    )

# Let the central slice theorem work its magic!
# Interpolate the 2D Fourier space grid from the transformed sinogram rows
fft2=scipy.interpolate.griddata(
    (srcy,srcx),
    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(np.real(fft2),vmin=-V,vmax=V)
plt.subplot(222)
plt.title("Recon (imag)")
plt.imshow(np.imag(fft2),vmin=-V,vmax=V)

# Show 2D FFT of target, just for comparison
expected_fft2=scipy.fftpack.fftshift(
    scipy.fftpack.fft2(
        scipy.fftpack.ifftshift(
            target
            )
        )
    )

plt.subplot(223)
plt.title("Expected (real)")
plt.imshow(np.real(expected_fft2),vmin=-V,vmax=V)
plt.subplot(224)
plt.title("Expected (imag)")
plt.imshow(np.imag(expected_fft2),vmin=-V,vmax=V)

# Transform from 2D Fourier space back to a reconstruction of the target
recon=np.real(
    scipy.fftpack.fftshift(
        scipy.fftpack.ifft2(
            scipy.fftpack.ifftshift(fft2)
            )
        )
    )

plt.figure()
plt.title("Reconstruction")
plt.imshow(recon,vmin=0.0,vmax=1.0)
plt.gray()

plt.show()

更新2013-02-17：如果您有足够的兴趣来遍历这堆东西，可以从这张海报的形式中找到更多自学程序的输出，这是自学程序的一部分。该存储库中的代码主体也可能很有趣（尽管请注意，代码并没有像上面的代码那样精简）。我可能会在某个时候尝试将其重新包装为IPython“笔记本”。

#2 楼

我不确切知道问题出在哪里，但是切片定理意味着这两种特殊情况应该是正确的：

找到这些停止点相等的点，从正弦图中向前移动，然后从生成的2D FFT中向下移动。

这看起来不正确：

fft2(target)[0] = fft(sinogram[270])
fft2(target)[:,0] = fft(sinogram[0])

您生成的2D FFT应该旋转90度？

我建议使用幅度和相位而不是实数和虚数，这样您可以更轻松地看到发生了什么事：

（白角是log(abs(0))的-inf，这不是问题）

#3 楼

我认为第一个解决方案无效的实际理论原因是因为旋转是相对于图像中心进行的，导致[S/2, S/2]的偏移，这意味着sinogram的每一行都没有从0到S，而是从-S/2到S/2。在您的示例中，偏移量实际上是offset = np.floor(S/2.)。请注意，这适用于S的偶数或奇数，并且等效于您在代码S/2中所做的操作（尽管更明确地避免了问题，例如S是float）。

我的猜测是，这种转变在傅立叶变换（FT）中引入的相位延迟是您在第二条消息中所讨论的内容：相位混乱，需要补偿为了能够应用Radon变换的反演。一个人需要深入研究这一理论，以便确定反函数按预期工作的确切条件。

为了补偿该偏移量，您可以像以前一样使用fftshift（将每行的中心放在开头，并且由于使用DFT实际上对应于计算S周期信号的傅立叶变换），则可以得到正确的东西），或者在计算sinogram FT时在复杂的傅立叶变换中显式地补偿此效果。在实践中，您可以删除ifftshift而不是：

sinogram_fft_rows=scipy.fftpack.fftshift(
    scipy.fftpack.fft(
        scipy.fftpack.ifftshift(
            sinogram,
            axes=1
            )
        ),
    axes=1
    )

，并用校正矢量乘以每一行：
这来自傅立叶变换属性，在考虑时移时（请查看FT维基百科页面上的“移位定理”，并申请等于- offset的移位-因为我们将图像放回了中心位置）。

同样，您可以对重建应用相同的策略，并通过在两个维度上但在另一个方向上的相位校正来替换fftshift（向后补偿）：

offset = np.floor(S/2.)
sinogram_fft_rows = scipy.fftpack.fftshift(
    scipy.fftpack.fft(sinogram, axis=1)
    * (np.exp(1j * 2.* np.pi * np.arange(S) * offset / S)),
    axes=1)

，这并不能改善您的解决方案，但可以从另一个角度阐明问题的理论。希望对您有所帮助！

此外，我不太喜欢使用fftshift，因为它会弄乱fft的计算方式。但是，在这种情况下，您需要在插值之前将FT的中心放在图像的中心，以获取fft2（或至少在设置r时要小心-这样就可以完全使其完全免费！）， fftshift真的很方便。但是，我更喜欢将该功能用于可视化目的，而不是在计算“核心”中使用。 :-)

最好的问候，

Jean-Louis

P.S .：您是否尝试过在不裁剪圆圈的情况下重建图像？可以在角落产生非常酷的模糊效果，在Instagram之类的程序中具有这样的功能会很好，不是吗？

编程黑洞网

用傅立叶方法进行层析成像重建的代码有什么问题？

评论

#1 楼

#2 楼

#3 楼