Radon变换原理及其应用

发表于 2023-12-28 更新于 2024-01-02 分类于 Image Processing ， Transformation

Radon变换是一种积分变换，用于将二维函数不同的投影的信号进行处理，并以此估计二维函数的信息。Radon变换在CT成像、地震信号处理等任务中有广泛的应用。

Radon变换的基本概念

Radon变换（中译名：拉东变换、雷登变换）是一种积分变换，通常用于从一系列角度的投影得到的1D信号去估计2D的信号分布情况。这个操作的一个重要和常见的应用就是CT扫描重建截面密度，从而了解人体器官的情况。Radon变换的任务设定如下图所示：

对于一个定义在2D空间的紧支撑函数\(f(x,y)\)，其每个坐标的取值代表该点的密度，我们希望得到这个二维密度函数在各个位置的密度，但是由于发射源和接收器无法以取点的方式查询2D空间，而只能将从source到detector的整条线上的密度进行积分。这个过程可以写为线积分的形式： \[ R_L = \int_L f(x, y) ds \] 通常我们有多个平行的source和detector之间的连线，因此一次会收集到多个detector的信息，从而得到一个1D的信号。这个1D信号对应着一个投影方向（即平行线的透射方向）。下面我们需要将线积分中的积分路径表示为直线，从而整理积分形式。

直线的黑塞标准式（Hesse normal form）

一般来说，我们用斜率-截距表达式来表示直线： \[ y = kx + b \] 但是这种表达方式对于垂直于（或者接近垂直于）\(x\)轴的直线来说，\(k\)值趋于无穷大，因此作为参数空间来说无界（unbound），为了解决这个问题，考虑如下的方式表示直线，该方式称为黑塞标准式（Hesse normal form）： \[ \vec{r} \cdot \vec{n_0} - d = 0 \] 其中，\(d\)表示原点到直线的距离，\(\vec{n_0}\)表示直线的单位法向量（可以看出，实际上这个形式可以推广到更高维度，即超平面的表示）。对于直线来说，单位法向量可用直线与\(x\)轴夹角\(\alpha\)的三角函数表示（根据几何关系很容易观察到）： \[ \vec{n_0} = (\cos\alpha, \sin\alpha) \]

\[ x\cos\alpha + y\sin\alpha = d \]

Radon变换的数学原理

将直线带入到线积分中，并且利用\(\delta\)函数，可以将线积分改写为： \[ R(\alpha, d) = \int\int f(x, y) \delta(d - x\cos\alpha - y\sin\alpha) dxdy \] 从参数\((x, y)\)到参数\((\alpha, d)\)的映射，就是Radon变换，它将2D空间中的坐标点及其信息转换到了角度和距离组成的2D空间。通常来说，我们会先固定一个\(\alpha\)，然后给一系列的\(d\)，从而得到以\(d\)为坐标轴的1D信号（如果觉得不好理解，可以想象一下做CT的过程的一个截面，先固定角度，然后所有的射线源一起发射，接收器接受到一个1D信号，然后角度旋转，继续上述过程...）。我们需要做的就是将这个以\(d\)为变量的1D信号与原始的2D信号\(f(x,y)\)建立联系。

实际上，对于上述线积分，等号两边同时取傅里叶变换（固定\(\alpha\)，对\(d\)取傅里叶变换），经过推导后可以得到： \[ F[f](\omega\cos\alpha_0,\omega\sin\alpha_0) = F[R(\alpha_0,d)] \] 左边为原2D函数\(f\)的二维傅里叶变换，根据自变量的2D坐标发现，这里等式左边规定了过零点且角度为\(\alpha_0\)的频谱上的直线各点的值，右边就是该角度下的1D积分结果的傅里叶变换。这样两者就建立了关系。

这个特性提示我们，可以反过来操作：先对不同的角度投影，计算积分结果的1D傅里叶变换，然后将该频谱结果放到对应角度上的二维频谱图上，当采样角度足够多，我们就可以得到原图的2D傅里叶变换。再通过傅里叶反变换就可重建原图！这就是Radon变换重建的整体流程思路。如下图所示：

Fourier 切片定理

傅里叶切片定理（Fourier slicing theorem），又称投影切片定理（ projection-slice theorem）或中心切片定理（central slice theorem），更简洁地说明了这个问题： \[ F_{1D}\cdot Proj_{1D} = Slice_{2D} \cdot F_{2D} \] 简单来说就是：对2D的傅里叶域取一个角度的切片，等价于先沿着该角度投影再做1D傅里叶变换。

实验示例

通过一个简单实验演示Radon变换的效果。

import numpy as np
import matplotlib.pyplot as plt

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

# 测试图为Shepp-Logan幻影图：https://en.wikipedia.org/wiki/Shepp-Logan_phantom
image = shepp_logan_phantom()
image = rescale(image, scale=0.4, mode='reflect', channel_axis=None)
# 测试不同的角度采样比例
for ratio in [0.1, 0.3, 0.5, 1.0]:
    theta = np.linspace(0., 180., int(ratio * max(image.shape)), endpoint=False)
    # radon transform
    sinogram = radon(image, theta=theta)
    # inverse radon transform
    recon = iradon(sinogram, theta=theta, filter_name='ramp')

    ax = plt.figure()
    ax.add_subplot(131)
    plt.imshow(image, aspect="auto")
    plt.title('image')
    ax.add_subplot(132)
    plt.imshow(sinogram, aspect="auto")
    plt.title(f'sinogram {ratio}')
    ax.add_subplot(133)
    plt.imshow(recon, aspect="auto")
    plt.title('reconstruction')

plt.show()