内插复数矩阵

Question

我想旋转一行复数（实际上是 Radon 变换的一行的一维 FFT），我在 Matlab 中使用 imrotate 但我不认为插值在做它应该是什么。

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

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

我需要获取 Radon 变换的每一行并根据其角度旋转它并将其放在二维矩阵中的相应角度。完成此操作后，我可以获得 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

旋转二维离散图像中的单条线真的很难。你总是得到一个粗略的近似值，插值并没有多大帮助。

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

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

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

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

旋转和 IFFT 运算也相互交换，因此以上内容等同于：

对于 Radon 变换中的每个投影：
- 应用 FFT
- 应用楔形过滤器
- 通过二维复杂图像的原点将其写成一行
- 应用二维 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行）的预期，直接在原点写入该行。生成的结果正确地再现了幻影。

内插复数矩阵

Interpolate matrix of complex numbers

python

image-processing

fft

interpolation

complex-numbers