前几天我学习了Radon变换并用Python做了一个简单的程序(见上一篇博文)，昨天看了一下逆Radon变换，尝试了一下简单的实现。我们可以通过对Sinogram图使用逆Radon变换来还原原始图像，最简单的方法是直接反投影，另一种方法是滤波反投影(这个方法是大部分商用CT机器使用的算法)。下面我介绍一下这两种方法，然后给出一个简单的Python实现：

直接反投影

直接反投影就是直接将投影值均匀回抹，然后将不同角度投影的回抹值相叠加得到原始的图片，代码如下：

import numpy as np

from scipy import ndimage

import matplotlib.pyplot as plt

import imageio

from cv2 import cv2

def IRandonTransform(image, steps):

#定义用于存储重建后的图像的数组

channels = len(image[0])

origin = np.zeros((steps, channels, channels))

for i in range(steps):

projectionValue = image[:, i]

projectionValueExpandDim = np.expand_dims(projectionValue, 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("radonshepplogan.png", 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()

注意：因为旋转角都是离散的，常常需要使用两个相邻角度的扫描数据来对中间的数据进行插值，但是上面的代码没有进行插值，而是直接进行了反投影。

跟上篇博文类似，我验证了两种情形，一种是Shepp-Logan模型，一种是直线情形，结果如下：

使用Shepp—Logan模型的投影值进行直接反投影得到的结果。

使用黑底白直线投影值进行直接反投影得到的结果。

以上就是直接回抹得到的结果，可以看出明显的伪影。因此通常需要对投影值进行滤波后再进行反投影。下面简单地介绍一下滤波反投影。

滤波反投影

顾名思义，滤波反投影就是对先对投影值进行滤波，然后利用的得到的值进行反投影，简单来说滤波就是使边缘(对应于图像的高频部分)更加明显。理论上来说，滤波函数应该是：

h(ω)=∣ω∣h(\omega) = |\omega|h(ω)=∣ω∣

但是这个是一个理想滤波器，没办法实现，因此需要考虑其他能够实现并且能够使得到的图像更加接近原始图像的滤波器。这里我仅介绍两种R—L滤波器和S—L滤波器，下面是这两种滤波器的具体形式：

hRL(nδ)={−14δ2,n=0,0,n为偶数,−1(nπδ)2,n为奇数.h_{RL}(n\delta) =

\begin{cases}

-\frac{1}{4\delta^{2}}, & \text{n=0,}\\

0, & \text{n为偶数,}\\

-\frac{1}{\left(n\pi\delta\right)^{2}}, & \text{n为奇数.}

\end{cases}hRL​(nδ)=⎩⎪⎨⎪⎧​−4δ21​,0,−(nπδ)21​,​n=0,n为偶数,n为奇数.​

hSL(nδ)=1π2δ2(4n2−1)h_{SL}(n\delta)=\frac{1}{\pi^{2}\delta^{2}\left(4n^{2}-1\right)}hSL​(nδ)=π2δ2(4n2−1)1​

利用以上两个滤波器和投影值进行卷积，然后再进行反投影，就可以得到得到原始的图像，大致上来说，就是在上面的代码中增加滤波器的构建和投影与滤波器的卷积过程，具体的代码如下：

import numpy as np

from scipy import ndimage

from scipy.signal import convolve

import matplotlib.pyplot as plt

import imageio

from cv2 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("radonshepplogan.png", 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()

下面我们来看一结果，仍旧是分别使用Shepp—Logan模型和黑底白直线进行测试：

对Shepp—Logan模型的投影值使用RL滤波器进行滤波反投影得到的结果。

对Shepp—Logan模型的投影值使用SL滤波器进行滤波反投影得到的结果。

对黑底白直线得到的投影值使用RL滤波器进行滤波反投影得到的结果。

对黑底白直线得到的投影值使用RL滤波器进行滤波反投影得到的结果。

通过对比直接反投影和滤波反投影，我们可以看到，因为高通滤波器的存在，重建的图像中原始图像的高频部分更加突出，边缘更加明显。需要注意的是，在滤波反投影过程中，我没有对投影信息进行插值，如果增加插值过程，应该能够获得更好的结果。希望大家多多批评指正！