拉东变换及其应用

• 1 算法背景
• 2 算法原理
• • 2.1 拉东变换
• 3 应用
• • 3.1 逆拉东变换
  • • 3.1.1 中心切片定理
    • 3.1.2 直接反投影
    • 3.1.3 滤波反投影（FBP）
• 4 测试代码
• • 4.1 使用skimage自带的接口
  • 4.2 使用理论编写

1 算法背景

拉东变换是由奥地利数学家约翰·拉东于1917年提出，目前被广泛的应用在断层扫描，其反变换可以从断层扫描的剖面图重建出投影前的函数。在数学上，拉东变换本质是一种积分变换，这个变换将二维平面函数 f 变换成一个定义在二维空间上的一个线性函数 Rf R_fRf​ 。 Rf R_fRf​的值为对函数$f$沿着直线$A$做积分的值，以下图为例：

2 算法原理

2.1 拉东变换

令密度函数 μ = μ ( x1, x2)μ=μ(x_1, x_2)μ=μ(x1​,x2​)，其定义域为 R2 R^2R2，令$R$为拉东变换的算子，则 Ɍ L μ ( x 1 , x 2 ) Ɍ_{Lμ(x1,x2)}ɌLμ(x1,x2)​表示沿着直线L对$\mu$的积分：
R Lμ = ∫ Lμ( x 1, x 2)d x 1d x 2R_{L\mu}=\int\limits_{L}^{} \mu(x_1,x_2)dx_1dx_2 RLμ​=L∫​μ(x1​,x2​)dx1​dx2​
一般会将直线通过法线式来表示： x1c o s θ + x2s i n θ − s = 0x_1cos\theta+x_2sin\theta-s=0x1​cosθ+x2​sinθ−s=0

故拉东变换可以采用下面的表示方式：
R Lμ =Rμ(s,θ)=∬μ( x 1, x 2)δ(s− x 1cosθ− x 2sinθ)d x 1d x 2 s∈(−∞,∞),θ∈(0,π)R_{L\mu}=R\mu(s, \theta)=\iint \mu(x_1,x_2)\delta(s-x_1cos\theta-x_2sin\theta)dx_1dx_2 \\s\in (-\infty, \infty),\theta\in(0, \pi) RLμ​=Rμ(s,θ)=∬μ(x1​,x2​)δ(s−x1​cosθ−x2​sinθ)dx1​dx2​s∈(−∞,∞),θ∈(0,π)

3 应用

CT成像：在医学图像成像时，是通过x射线穿过人体来实现的，射线穿过人体后会衰减，然后被仪器测量到，也就是说我们已经知道了拉东变换的结果，需要通过这个结果还原出射线穿过的人体剖面，即拉东逆变换。

3.1 逆拉东变换

3.1.1 中心切片定理

中心切片定理为拉东逆变换提供了数学上证明：可以根据投影值完全可以重建原始图像。该定理可以表示为密度函数 μ( x 1, x 2)\mu(x_1, x_2) μ(x1​,x2​) 沿某个方向的投影函数 $\rho (s,\theta )$ 的傅里叶变换等于函数 μ( x 1, x 2)\mu(x_1, x_2) μ(x1​,x2​)的二维傅里叶变换沿探测器平行方向过原点的片段，如下图：

证明如下，在角度$\theta$已知的情况下对$p(s,\theta )$进行傅里叶变换:F(p(s,θ))= ∫ −∞∞p(s,θ) e −jws dsF(p(s,\theta))=\int_{-\infty}^{\infty}p(s,\theta)e^{-jws}ds F(p(s,θ))=∫−∞∞​p(s,θ)e−jwsds
在已知 p ( s , θ ) = ∬ μ ( x1, x2) δ ( s − x1c o s θ − x2s i n θ ) d x1d x2 p(s, \theta)=\iint \mu(x_1,x_2)\delta(s-x_1cos\theta-x_2sin\theta)dx_1dx_2p(s,θ)=∬μ(x1​,x2​)δ(s−x1​cosθ−x2​sinθ)dx1​dx2​的情况下，可得：
F(p(s,θ))= ∫ −∞∞p(s,θ) e −jws ds= ∫ −∞∞(∬μ( x 1, x 2)δ(s− x 1cosθ− x 2sinθ)d x 1d x 2) e −jws ds= ∬μ( x 1, x 2)[δ(s− x 1cosθ− x 2sinθ) e −jws ds]d x 1d x 2= ∬μ( x 1, x 2) e −jw( x 1cosθ+ x 2sinθ) d x 1d x 2F(p(s,\theta))=\int_{-\infty}^{\infty}p(s,\theta)e^{-jws}ds=\int_{-\infty}^{\infty}(\iint \mu(x_1,x_2)\delta(s-x_1cos\theta-x_2sin\theta)dx_1dx_2)e^{-jws}ds=\\\iint \mu(x_1,x_2)[\delta(s-x_1cos\theta-x_2sin\theta)e^{-jws}ds]dx_1dx_2=\\\iint \mu(x_1,x_2)e^{-jw(x_1cos\theta+x_2sin\theta)}dx_1dx_2 F(p(s,θ))=∫−∞∞​p(s,θ)e−jwsds=∫−∞∞​(∬μ(x1​,x2​)δ(s−x1​cosθ−x2​sinθ)dx1​dx2​)e−jwsds=∬μ(x1​,x2​)[δ(s−x1​cosθ−x2​sinθ)e−jwsds]dx1​dx2​=∬μ(x1​,x2​)e−jw(x1​cosθ+x2​sinθ)dx1​dx2​
中心切片定理一般不用于投影重建。通过下图可知：
当往 wx− wy w_x-w_ywx​−wy​平面一条一条地添加“中心片段”时，中心片段在 wx− wy w_x-w_ywx​−wy​平面的原点密度高于远离原点区域的密度。这就意味着高频区域的大片面积是没有中心片段覆盖的，但这些没被覆盖的地方只能通过插值数据。但凡是进行插值，一定会引入大量误差，因此用中心切片定理重建图像时，往往图像的质量不好。

3.1.2 直接反投影

直接反投影即将投影值均匀回抹，然后将不同角度的回抹值叠加得到重建图像。


反投影图是离散叠加的，显然在中心处信号集中，边缘处信号稀疏，因此在最后需要在空缺的地方进行插值，才能得到最终的原图像。

3.1.3 滤波反投影（FBP）

当前主流的反投影重建算法主要有滤波反投（FBP）和迭代重建算法。滤波反投影的主要缺点是噪声大，信噪比低，但由于处理数据较少所以重建速度快。迭代重建主要是计算量特别大，重建速度慢，但是图像信噪比较高，图像质量相对较好。
滤波反投影算法的步骤：

1. 计算每一个投影的一维傅里叶变换；

2. 用滤波函数乘以每一个傅里叶变换；

3. 得到每一个滤波后的变换的一维反傅里叶变换；

4. 对步骤3得到的所以一维反变换积分（求和）
   可以通过在频域对投影函数乘上一个高通滤波器$|w|$的方式来抵消直接反投影算法带来的误差。但这是一个理想滤波器，没办法实现，因此需要考虑其他能够实现并且能够使得到的图像更加接近原始图像的滤波器。
   R-L滤波器是使用窗函数对斜坡滤波器进行截断产生的。如下图所示：
   
   离散化的函数表达式如下：
   

下面为采用R-L滤波器的结果：

4 测试代码

4.1 使用skimage自带的接口

from cgitb import greyimport cv2 as cvimport numpy as npimport matplotlib.pyplot as pltfrom skimage import transformfrom skimage.transform import radonfrom skimage import ioimage = io.imread('test.jpg', as_gray=grey)image = transform.resize(image, (image.shape[0], image.shape[0]))##将图像转换为正方形，方便后续滤波计算# 图像坐标转换为(theta,p)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., 180., max(image.shape), endpoint=False)sinogram = radon(image, theta=theta)print(np.shape(sinogram))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()cv.imwrite("ladon_sinogram.jpg", sinogram)#反变化结果from skimage.transform import iradonsize = np.shape(image)reconstruction_fbp = iradon(sinogram, theta=theta, filter_name='cosine')error = reconstruction_fbp - imageprint(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()

4.2 使用理论编写

import numpy as npfrom scipy import ndimagefrom scipy.signal import convolveimport matplotlib.pyplot as pltimport imageioimport cv2#两种滤波器的实现def RLFilter(N, d):filterRL = np.zeros((N,))for i in range(N):filterRL[i] = - 1.0 / np.power((i - N / 2) * np.pi * d, 2.0)if np.mod(i - N / 2, 2) == 0:filterRL[i] = 0filterRL[int(N/2)] = 1 / (4 * np.power(d, 2.0))return filterRLdef SLFilter(N, d):filterSL = np.zeros((N,))for i in range(N):#filterSL[i] = - 2 / (np.power(np.pi, 2.0) * np.power(d, 2.0) * (np.power((4 * (i - N / 2)), 2.0) - 1))filterSL[i] = - 2 / (np.pi**2.0 * d**2.0 * (4 * (i - N / 2)**2.0 - 1))return filterSLdef IRandonTransform(image, steps):#定义用于存储重建后的图像的数组channels = len(image[0])origin = np.zeros((steps, channels, channels))filter = RLFilter(channels, 1)# filter = SLFilter(channels, 1)for i in range(steps):projectionValue = image[:, i]projectionValueFiltered = convolve(filter, projectionValue, "same")projectionValueExpandDim = np.expand_dims(projectionValueFiltered, axis=0)projectionValueRepat = projectionValueExpandDim.repeat(channels, axis=0)origin[i] = ndimage.rotate(projectionValueRepat, i*180/steps, reshape=False).astype(np.float64)iradon = np.sum(origin, axis=0)return iradon#读取图片#image = cv2.imread('straightLine.png', cv2.IMREAD_GRAYSCALE)image = cv2.imread("radon.jpg", cv2.IMREAD_GRAYSCALE)iradon = IRandonTransform(image, len(image[0]))#绘制原始图像和对应的sinogram图plt.subplot(1, 2, 1)plt.imshow(image, cmap='gray')plt.subplot(1, 2, 2)plt.imshow(iradon, cmap='gray')plt.show()

测试图像：