拉东变换介绍和代码实现

一 、正向变换

作为我们的原始图像，我们将使用 Shepp-Logan 幻影。在计算拉东变换时，我们需要确定要使用多少个投影角度。经验法则表明，投影的数量应该与对象横跨的像素数量大致相同（要了解原因，请考虑在重建过程中必须确定多少个未知像素值，并将此值与投影提供的测量值进行比较），我们在这里遵循此规则。以下是原始图像及其拉东变换，通常称为其正弦图

import numpy as np
import matplotlib.pyplot as plt
from skimage.data import shepp_logan_phantom
from skimage.transform import radon, rescale

# 导入图像处理库skimage中的数据集和变换函数
from skimage import data

# 生成并缩放Shepp-Logan伪影图像
image = data.shepp_logan_phantom()
image = rescale(image, scale=0.4, mode='reflect', channel_axis=None)

# 创建子图用于显示原始图像和Radon变换结果
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5))

# 设置第一个子图标题并显示原始图像
ax1.set_title("Original")
ax1.imshow(image, cmap=plt.cm.Greys_r, interpolation='nearest')

# 定义Radon变换的投影角度,image.shape即二维图像的宽度和高度
theta = np.linspace(0.0, 180.0, max(image.shape), endpoint=False)
# 执行Radon变换并生成投影图
sinogram = radon(image, theta=theta)

# 设置第二个子图标题并显示Radon变换结果
ax2.set_title("Radon transform\n(Sinogram)")
ax2.set_xlabel("Projection angle (deg)")
ax2.set_ylabel("Projection position (pixels)")
# 显示投影图并调整显示范围和纵横比
ax2.imshow(
    sinogram,
    cmap=plt.cm.Greys_r,
    extent=(-np.pi, 180.0, -1, sinogram.shape[0]),
    aspect='auto',
)

# 调整子图布局以适应显示区域
fig.tight_layout()
# 显示图像
plt.show()

 输出结果：

二、使用滤波反投影 (FBP) 进行重建

波反投影的数学基础是傅里叶切片定理[2]。它使用投影的傅里叶变换和傅里叶空间中的插值来获得图像的 2D 傅里叶变换，然后将其反转以形成重建图像。滤波反投影是执行逆拉东变换的最快速方法之一。FBP 的唯一可调参数是应用于傅里叶变换投影的滤波器。它可用于抑制重建中的高频噪声。skimage 提供了“ramp”、“shepp-logan”、“cosine”、“hamming”和“hann”滤波器

# 导入matplotlib.pyplot库用于绘图
# 导入skimage.transform.radon_transform模块中的_get_fourier_filter函数
import matplotlib.pyplot as plt
from skimage.transform.radon_transform import _get_fourier_filter

# 定义一个包含不同滤波器名称的列表
filters = ['ramp', 'shepp-logan', 'cosine', 'hamming', 'hann']

# 遍历滤波器列表
for ix, f in enumerate(filters):
    # 调用_get_fourier_filter函数获取指定滤波器的响应
    response = _get_fourier_filter(2000, f)
    # 使用plt.plot绘制滤波器响应的频率图，并添加图例
    plt.plot(response, label=f)

# 设置x轴的显示范围为0到1000
plt.xlim([0, 1000])
# 设置x轴的标签为频率
plt.xlabel('frequency')
# 显示图例
plt.legend()
# 显示绘制的图形
plt.show()

输出结果：

三、 应用“ramp”滤波器执行逆拉东变换

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

# 读取并预处理图像
image = shepp_logan_phantom()
image = rescale(image, scale=0.4, mode='reflect', channel_axis=None)  # 对图像进行缩放和反射处理

# 生成投影角度和投影图像
theta = np.linspace(0.0, 180.0, max(image.shape), endpoint=False)
sinogram = radon(image, theta=theta)  # 使用Radon变换生成投影图像

# 定义滤波器列表
filters = ['ramp', 'shepp-logan', 'cosine', 'hamming', 'hann']  # 定义不同的滤波器名称

# 使用iradon函数进行滤波反投影重建
reconstruction_fbp = iradon(sinogram, theta=theta, filter_name='ramp')  # 使用滤波反投影重建图像
error = reconstruction_fbp - image  # 计算重建图像与原始图像之间的误差
print(f'FBP rms reconstruction error: {np.sqrt(np.mean(error**2)):.3g}')  # 打印重建误差的均方根值

# 设置图像显示参数
imkwargs = dict(vmin=-0.2, vmax=0.2)  # 设置图像显示的最小和最大值

# 创建画布和子图
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5), sharex=True, sharey=True)  # 创建包含两个子图的画布

# 显示重建图像和重建误差图像
ax1.imshow(reconstruction_fbp, cmap=plt.cm.Greys_r)  # 显示重建图像
ax2.imshow(error, cmap=plt.cm.Greys_r, **imkwargs)  # 显示重建误差图像

# 显示图像
plt.show()  # 显示画布和子图中的图像

生成结果如下图：