拉东变换及其应用

大家好，欢迎来到IT知识分享网。

1 算法背景

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

2 算法原理

2.1 拉东变换

令密度函数$\mu =\mu (x_1,x_2)$，其定义域为$R^2$，令$R$为拉东变换的算子，则${Ɍ}_{L\mu (x1,x2)}$表示沿着直线L对$\mu$的积分：
$$R_{L\mu }=\underset{L}{\int }\mu (x_1,x_2)d{x}_{1}d{x}_2$$
一般会将直线通过法线式来表示：$x_{1}cos\theta +x_{2}sin\theta -s=0$

3 应用

3.1 逆拉东变换

3.1.1 中心切片定理

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

证明如下，在角度$\theta$已知的情况下对$p(s,\theta )$进行傅里叶变换:$$F(p(s,\theta ))={\int }_{-\infty }^{\infty }p(s,\theta )e^{-jws}ds$$
在已知$p(s,\theta )=\iint \mu (x_1,x_2)\delta (s-x_{1}cos\theta -x_{2}sin\theta )d{x}_{1}d{x}_2$的情况下，可得：
$$\begin{gathered} F(p(s,\theta ))={\int }_{-\infty }^{\infty }p(s,\theta )e^{-jws}ds={\int }_{-\infty }^{\infty }(\iint \mu (x_1,x_2)\delta (s-x_{1}cos\theta -x_{2}sin\theta )d{x}_{1}d{x}_2)e^{-jws}ds= \\ \iint \mu (x_1,x_2)[\delta (s-x_{1}cos\theta -x_{2}sin\theta )e^{-jws}ds]d{x}_{1}d{x}_2= \\ \iint \mu (x_1,x_2)e^{-jw(x_{1}cos\theta +x_{2}sin\theta )}d{x}_{1}d{x}_2 \end{gathered}$$
中心切片定理一般不用于投影重建。通过下图可知：
当往$w_x-w_y$平面一条一条地添加“中心片段”时，中心片段在$w_x-w_y$平面的原点密度高于远离原点区域的密度。这就意味着高频区域的大片面积是没有中心片段覆盖的，但这些没被覆盖的地方只能通过插值数据。但凡是进行插值，一定会引入大量误差，因此用中心切片定理重建图像时，往往图像的质量不好。

3.1.2 直接反投影

3.1.3 滤波反投影（FBP）

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

4 测试代码

4.1 使用skimage自带的接口

from cgitb import grey import cv2 as cv import numpy as np import matplotlib.pyplot as plt from skimage import transform from skimage.transform import radon from skimage import io image = 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 iradon size = np.shape(image) reconstruction_fbp = iradon(sinogram, theta=theta, filter_name='cosine') error = reconstruction_fbp - image print(f'FBP rms reconstruction error: { 
     np.sqrt(np.mean(error2)):.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 np from scipy import ndimage from scipy.signal import convolve import matplotlib.pyplot as plt import imageio import 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] = 0 filterRL[int(N/2)] = 1 / (4 * np.power(d, 2.0)) return filterRL def 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.pi2.0 * d2.0 * (4 * (i - N / 2)2.0 - 1)) return filterSL def 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()