备注

转到结尾 以下载完整示例代码。或在浏览器中通过 Binder 运行此示例。

Radon 变换#

在计算机断层扫描中，断层重建问题是从一组投影中获得断层切片图像 [1]。投影是通过在感兴趣的2D物体上绘制一组平行光线形成的，将物体对比度沿每条光线的积分分配给投影中的单个像素。2D物体的单个投影是一维的。为了实现物体的计算机断层扫描重建，必须获取多个投影，每个投影对应于光线相对于物体的不同角度。在多个角度收集的投影称为正弦图，它是原始图像的线性变换。

逆Radon变换用于计算机断层扫描中，从测量的投影（正弦图）重建二维图像。逆Radon变换的实际精确实现并不存在，但有几种良好的近似算法可用。

由于逆拉东变换从一组投影中重建物体，因此（正向）拉东变换可以用来模拟断层扫描实验。

此脚本执行拉东变换以模拟断层扫描实验，并基于模拟生成的正弦图重建输入图像。比较了两种执行逆拉东变换和重建原始图像的方法：滤波反投影（FBP）和同时代数重建技术（SART）。

关于断层重建的更多信息，请参见：

正向变换#

作为我们的原始图像，我们将使用 Shepp-Logan 模型。在计算 Radon 变换时，我们需要决定使用多少个投影角度。根据经验，投影的数量应与物体横跨的像素数量大致相同（要了解为什么如此，请考虑重建过程中必须确定的未知像素值的数量，并将其与投影提供的测量数量进行比较），我们在这里遵循这一规则。以下是原始图像及其 Radon 变换，通常称为其 正弦图：

import numpy as np
import matplotlib.pyplot as plt

from skimage.data import shepp_logan_phantom
from skimage.transform import radon, rescale

image = shepp_logan_phantom()
image = rescale(image, scale=0.4, mode='reflect', channel_axis=None)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5))

ax1.set_title("Original")
ax1.imshow(image, cmap=plt.cm.Greys_r)

theta = np.linspace(0.0, 180.0, max(image.shape), endpoint=False)
sinogram = radon(image, theta=theta)
dx, dy = 0.5 * 180.0 / max(image.shape), 0.5 / sinogram.shape[0]
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=(-dx, 180.0 + dx, -dy, sinogram.shape[0] + dy),
    aspect='auto',
)

fig.tight_layout()
plt.show()
Original, Radon transform (Sinogram)

使用滤波反投影（FBP）进行重建#

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

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):
    response = _get_fourier_filter(2000, f)
    plt.plot(response, label=f)

plt.xlim([0, 1000])
plt.xlabel('frequency')
plt.legend()
plt.show()
plot radon transform

应用带有’ramp’滤波器的逆拉东变换，我们得到：

from skimage.transform import 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.set_title("Reconstruction\nFiltered back projection")
ax1.imshow(reconstruction_fbp, cmap=plt.cm.Greys_r)
ax2.set_title("Reconstruction error\nFiltered back projection")
ax2.imshow(reconstruction_fbp - image, cmap=plt.cm.Greys_r, **imkwargs)
plt.show()
Reconstruction Filtered back projection, Reconstruction error Filtered back projection
FBP rms reconstruction error: 0.0283

使用同步代数重建技术进行重建#

代数重建技术用于断层扫描基于一个简单的想法：对于一个像素化的图像，特定投影中单条射线的值是射线穿过物体时经过的所有像素值的总和。这是正向拉东变换的一种表达方式。逆拉东变换可以被表述为一个（大）线性方程组。由于每条射线只穿过图像中一小部分像素，这个方程组是稀疏的，允许使用稀疏线性系统的迭代求解器来处理这个方程组。一种迭代方法特别受欢迎，即Kaczmarz方法[3]_，它具有这样的特性：解将趋近于方程组的最小二乘解。

将重建问题公式化为一组线性方程与迭代求解器的结合，使得代数技术相对灵活，因此某些形式的先验知识可以相对容易地被纳入。

skimage 提供了代数重建技术中较为流行的一种变体：同步代数重建技术（SART）[4]_。它使用Kaczmarz方法作为迭代求解器。通常在一次迭代中就能获得良好的重建效果，使得该方法在计算上非常有效。运行一次或多次额外的迭代通常会改善锐利、高频特征的重建，并减少均方误差，但代价是增加了高频噪声（用户需要根据当前问题的具体情况决定最佳迭代次数。skimage 中的实现允许在重建过程中提供形式为重建值上下阈值的先验信息。

from skimage.transform import iradon_sart

reconstruction_sart = iradon_sart(sinogram, theta=theta)
error = reconstruction_sart - image
print(
    f'SART (1 iteration) rms reconstruction error: ' f'{np.sqrt(np.mean(error**2)):.3g}'
)

fig, axes = plt.subplots(2, 2, figsize=(8, 8.5), sharex=True, sharey=True)
ax = axes.ravel()

ax[0].set_title("Reconstruction\nSART")
ax[0].imshow(reconstruction_sart, cmap=plt.cm.Greys_r)

ax[1].set_title("Reconstruction error\nSART")
ax[1].imshow(reconstruction_sart - image, cmap=plt.cm.Greys_r, **imkwargs)

# Run a second iteration of SART by supplying the reconstruction
# from the first iteration as an initial estimate
reconstruction_sart2 = iradon_sart(sinogram, theta=theta, image=reconstruction_sart)
error = reconstruction_sart2 - image
print(
    f'SART (2 iterations) rms reconstruction error: '
    f'{np.sqrt(np.mean(error**2)):.3g}'
)

ax[2].set_title("Reconstruction\nSART, 2 iterations")
ax[2].imshow(reconstruction_sart2, cmap=plt.cm.Greys_r)

ax[3].set_title("Reconstruction error\nSART, 2 iterations")
ax[3].imshow(reconstruction_sart2 - image, cmap=plt.cm.Greys_r, **imkwargs)
plt.show()
Reconstruction SART, Reconstruction error SART, Reconstruction SART, 2 iterations, Reconstruction error SART, 2 iterations
SART (1 iteration) rms reconstruction error: 0.0329
SART (2 iterations) rms reconstruction error: 0.0214

脚本总运行时间： (0 分钟 0.681 秒)

由 Sphinx-Gallery 生成的图库

On this page
显示源代码