- 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 变换成一个定义在二维空间上的一个线性函数 。
的值为对函数
沿着直线
做积分的值，以下图为例：
2 算法原理
2.1 拉东变换
令密度函数，其定义域为
，令
为拉东变换的算子，则
表示沿着直线L对
的积分：
一般会将直线通过法线式来表示：
故拉东变换可以采用下面的表示方式：
3 应用
CT成像：在医学图像成像时，是通过x射线穿过人体来实现的，射线穿过人体后会衰减，然后被仪器测量到，也就是说我们已经知道了拉东变换的结果，需要通过这个结果还原出射线穿过的人体剖面，即拉东逆变换。
3.1 逆拉东变换
3.1.1 中心切片定理
中心切片定理为拉东逆变换提供了数学上证明：可以根据投影值完全可以重建原始图像。该定理可以表示为密度函数 沿某个方向的投影函数
的傅里叶变换等于函数
的二维傅里叶变换沿探测器平行方向过原点的片段，如下图：
3.1.2 直接反投影
直接反投影即将投影值均匀回抹，然后将不同角度的回抹值叠加得到重建图像。
3.1.3 滤波反投影（FBP）
当前主流的反投影重建算法主要有滤波反投（FBP）和迭代重建算法。滤波反投影的主要缺点是噪声大，信噪比低，但由于处理数据较少所以重建速度快。迭代重建主要是计算量特别大，重建速度慢，但是图像信噪比较高，图像质量相对较好。
滤波反投影算法的步骤：
-
计算每一个投影的一维傅里叶变换；
-
用滤波函数乘以每一个傅里叶变换；
-
得到每一个滤波后的变换的一维反傅里叶变换；
-
对步骤3得到的所以一维反变换积分（求和）
可以通过在频域对投影函数乘上一个高通滤波器的方式来抵消直接反投影算法带来的误差。但这是一个理想滤波器，没办法实现，因此需要考虑其他能够实现并且能够使得到的图像更加接近原始图像的滤波器。
R-L滤波器是使用窗函数对斜坡滤波器进行截断产生的。如下图所示：
下面为采用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(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 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.pi**2.0 * d**2.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()
测试图像：
文章出处登录后可见！
已经登录？立即刷新
