2020年8月18日补充 我后面的一篇博文实现了使用线性插值对投影值进行插值之后再进行滤波反投影的过程。

https://blog.csdn.net/hsyxxyg/article/details/108077396

插值之后得到的结果更好，大家可以去看那篇文章。

前几天我学习了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 cv2def IRandonTransform(image, steps):#定义用于存储重建后的图像的数组channels = len(image[0])origin = np.zeros((steps, channels, channels))for i in range(steps):#传入的图像中的每一列都对应于一个角度的投影值#这里用的图像是上篇博文里得到的Radon变换后的图像裁剪后得到的projectionValue = image[:, i]#这里利用维度扩展和重复投影值数组来模拟反向均匀回抹过程projectionValueExpandDim = np.expand_dims(projectionValue, axis=0)projectionValueRepeat = projectionValueExpandDim.repeat(channels, axis=0)origin[i] = ndimage.rotate(projectionValueRepeat, i*180/steps, reshape=False).astype(np.float64)#各个投影角度的投影值已经都保存在origin数组中，只需要将它们相加即可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] = 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)projectionValueRepeat = projectionValueExpandDim.repeat(channels, axis=0)origin[i] = ndimage.rotate(projectionValueRepeat, 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滤波器进行滤波反投影得到的结果。

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

测试代码时使用的图片

注意： 上面已经介绍了逆Radon变换的结果，下面是用到的原始图片（也就是上面的几张图片里左侧的图片），这些图片是我上篇博文里得到的结果，大家可以直接保存来测试上述代码或者其他的自己写的逆Radon变换代码。

Shepp—Logan模型的Radon变换的结果。

黑底白直线的Radon变换的结果。

Python实现逆Radon变换——直接反投影和滤波反投影相关推荐

pythonsl火车加字_荐Python实现Radon变换——直接反投影和滤波反投影
前几天我学习了Radon变换并用Python做了一个简单的程序(见上一篇博文),昨天看了一下逆Radon变换,尝试了一下简单的实现.我们可以通过对Sinogram图使用逆Radon变换来还原原始图像, ...
Python实现平行束滤波反投影——Inverse Radon Transformation
参考博文进行了平行束滤波反投影的修改,将时域滤波修改为频域滤波,重建后消除原博文中图像的竖条状伪影. https://blog.csdn.net/hsyxxyg/article/details/106 ...
【youcans 的 OpenCV 例程 200 篇】112. 滤波反投影重建图像
欢迎关注『youcans 的 OpenCV 例程 200 篇』系列,持续更新中欢迎关注『youcans 的 OpenCV学习课』系列,持续更新中 [youcans 的 OpenCV 例程 2 ...
【转】CT图像重构方法详解——傅里叶逆变换法、直接反投影法、滤波反投影法
转自:CT图像重构方法详解--傅里叶逆变换法.直接反投影法.滤波反投影法_Absolute Zero-CSDN博客_反投影法绪在做CT图像处理的时候遇到很多问题,对于滤波反变换有许多细 ...
CT图像重构方法详解——傅里叶逆变换法、直接反投影法、滤波反投影法
绪在做CT图像处理的时候遇到很多问题,对于滤波反变换有许多细节存在疑问,经过多天查找资料和利用MATLAB程序一步步实现后终于豁然开朗,于是想要总结成文,作为笔记方便今后查看.文中若有错误欢迎指出! ...
【OpenCV 例程 300 篇】112. 滤波反投影重建图像
专栏地址:『youcans 的 OpenCV 例程 300篇 - 总目录』 [第 7 章:图像复原与重建] 110. 投影和雷登变换 111. 雷登变换反投影重建图像 112. 滤波反投影重建图像 [ ...
运用滤波反投影的方法对图像进行重建matlab仿真
目录 1.算法描述 2.仿真效果预览 3.MATLAB部分代码预览 4.完整MATLAB程序 1.算法描述直接由正弦图得到反投影图像,会存在严重的模糊,这是早期 CT 系统所存在的问题.傅立叶中心切 ...
【CT算法，radon变换】基于MATLAB的CT算法，radon变换的三维建模仿真
1.软件版本 MATLAB2021a 2.本算法理论知识 1.输入:T(x,y,z) 使用stl读取函数完成T的导入工作 2.做Radon变换,得投影图:P 正常Radon变换即可. 3.对P:应用斜 ...
python画饼图_百度飞桨PaddlePaddle之[Python小白逆袭大神]7天训练营
第三次参加百度的7天训练营了这次参加的主题是[Python小白逆袭大神],不过你别看是小白逆势...除非你一开始参加就逆袭完,不然你真的是python小白,这个课程还是有难难度的. 说一下个训练营的 ...
反Radon变换 C++实现
文章目录 Radon变换原理代码实现对图像做滤波对图像做反投影创建新的dib对象并写入数据完整代码 Radon变换原理百度百科里面说的蛮清楚的,可以自己看一下. 频域添加了Ramp Fil ...

Python实现逆Radon变换——直接反投影和滤波反投影

直接反投影

滤波反投影

测试代码时使用的图片

Python实现逆Radon变换——直接反投影和滤波反投影相关推荐

最新文章

热门文章