插值复数矩阵

Question

我想旋转一行复数（实际上是 Radon 变换的一行的 1D FFT），我在 Matlab 中使用imrotate，但我认为插值没有做它应该做的事情。

目标是使用投影切片定理重现从氡到图像空间的转换。

^{（图片来自维基百科）}

我需要获取每一行的 Radon 变换并根据它的角度旋转它，并将它放在一个 2D 矩阵中的相应角度。完成此操作后，我可以获得 2D ifft2 来恢复图像（理论上）。这是目标。有人可以帮忙吗？

我想使用imrotate，但也许这不是正确的做法？目标是将 Radon 变换的 FFT 行映射到它们在圆圈中的正确位置，如上图所示。

这是旋转和最近邻插值的实际结果。右边的结果应该是通常的 SheppLogan 幻象。

import numpy as np
import matplotlib.pyplot as plt
from skimage.io import imread
from skimage.data import shepp_logan_phantom
from skimage.transform import radon, rescale
from skimage.transform import iradon
from skimage.transform import rotate
import cv2 


x=shepp_logan_phantom()
x=cv2.resize(x, (128,128), interpolation = cv2.INTER_AREA)
theta=np.linspace(0,180,len(x))
R=radon(x,theta)

temp_=np.zeros((128,128)).astype(np.complex128)
fullFft2D=np.zeros((128,128)).astype(np.complex128)

for i in range(len(theta)):
    temp_[63,:]=np.fft.fftshift(np.fft.fft(R[:,i])).T
    fft_real=rotate(np.real(temp_),theta[i],order=0)
    fft_imag=rotate(np.imag(temp_),theta[i],order=0)
    fullFft2D += fft_real+1j*fft_real
    temp_=np.zeros((128,128)).astype(np.complex128)


plt.imshow(np.fft.fftshift(np.abs(np.fft.ifft2(np.fft.ifftshift(fullFft2D)))))

我已经实现了你（@Luengo）所说的：

res=np.zeros((128,128))
tmp_=np.zeros((128,128)).astype(np.complex128)
for i in range(len(theta)):
    kspace_row = np.fft.fftshift(np.fft.fft(R[:,i])).T
    tmp_[63,:] = kspace_row
    res +=  rotate(np.abs(np.fft.ifft(np.fft.fftshift(tmp_))),-theta[i])

plt.imshow(res)

但它不起作用（我可能错过了什么？）

Answer 1

在 2D 离散图像中旋转单行非常困难。你总是得到一个粗略的近似值，插值没有多大帮助。

您打算遵循的过程是（我添加了过滤）：

• 对于 Radon 变换中的每个投影：
  • 应用 FFT
  • 应用楔形过滤器
  • 通过二维复杂图像的原点将其写成一行
  • 旋转此图像以匹配投影方向
  • 通过求和将结果累加到输出频域图像中
• 将2D IFFT应用于频域图像，得到重构图像

因为我们知道IFFT运算与求和对易，所以可以将IFFT运算移入循环：

• 对于 Radon 变换中的每个投影：
  • 应用 FFT
  • 应用楔形过滤器
  • 通过二维复杂图像的原点将其写成一行
  • 旋转此图像以匹配投影方向
  • 应用 2D IFFT
  • 通过求和将结果累加到输出的空间域图像中

旋转和 IFFT 操作也是通勤的，所以上面的内容等同于：

• 对于 Radon 变换中的每个投影：
  • 应用 FFT
  • 应用楔形过滤器
  • 通过二维复杂图像的原点将其写成一行
  • 应用 2D IFFT
  • 旋转此图像以匹配投影方向
  • 通过求和将结果累加到输出的空间域图像中

在后一种情况下，我们正在旋转一个平滑的空间域图像；它不是在否则为空的图像中绘制的一条线，它是一个完全带限制的函数，可以正确插值。在这种情况下，旋转结果要好得多。

后一个过程几乎与反投影算法所做的相同。我们可以进一步意识到，具有通过原点的单行数据（图像的其余部分全为零）的图像的 2D IFFT 与采用 1D IFFT 并在图像的所有行中复制它是相同的。这节省了相当多的计算量。

---

这是一些代码。第一种方法是（对 OP 的代码进行了一些修复，但输出仍然无法识别！）：

import numpy as np
import matplotlib.pyplot as plt
from skimage.io import imread
from skimage.data import shepp_logan_phantom
from skimage.transform import radon, rescale
from skimage.transform import iradon
from skimage.transform import rotate
import cv2 

x = shepp_logan_phantom()
x = cv2.resize(x, (128,128), interpolation = cv2.INTER_AREA)
theta = np.linspace(0, 180, len(x), endpoint=False)
R = radon(x, theta)

filter = np.abs(np.fft.fftfreq(128))

fullFft2D = np.zeros((128,128)).astype(np.complex128)
for i in range(len(theta)):
    temp_ = np.zeros((128,128)).astype(np.complex128)
    temp_[64,:] = np.fft.fftshift(filter * np.fft.fft(R[:,i]))
    fft_real = rotate(np.real(temp_), theta[i], order=0, center=(64,64))
    fft_imag = rotate(np.imag(temp_), theta[i], order=0, center=(64,64))
    fullFft2D += fft_real + 1j*fft_imag

y = np.fft.ifft2(np.fft.ifftshift(fullFft2D))
plt.imshow(np.fft.fftshift(y.real)); plt.show()

修复包括：(1) 大小为 128 的 fftshifted 频域中的原点为 64，而不是 63。(2) 明确地围绕原点执行旋转。 (3) OP 有一个错字：fft_real + 1j* fft_real。 (4) 增加了楔形滤波。 (5) 不包括 Radon 变换中的 180 度（因为它与 0 度相同）。 (6) 使用IFFT的实部，而不是绝对值。

在通过频域进行计算时，如果您期望得到实值结果，但得到非平凡（平凡==几乎为零）的虚部，则有问题。在上面的代码中，虚构的组件是不平凡的。这是无法正确插入的数据旋转的结果。轮换只是破坏了成功的变化。

后一种方法是：

y = np.zeros((128,128))
for i in range(len(theta)):
    tmp_ = np.zeros((128,128)).astype(np.complex128)
    tmp_[0,:] = filter * np.fft.fft(R[:,i])
    y += rotate(np.fft.ifft2(tmp_).real, -theta[i], center=(64,64))

plt.imshow(y); plt.show()

由于我们不需要使用fftshift，这段代码被稍微简化了，我们可以按照FFT（第0行）的预期直接在原点写入行。产生的结果正确地再现了幻影。