相关文章：编程记录——研究一下python对shepp_logan体模数据实现radon变换

参考博客：
CT典型数据——shepp_logan体模数据的生成 python版本
Python实现逆Radon变换——直接反投影和滤波反投影
主要是对上述第二个链接代码进行了少量改动总结。

一、补充

关于ndimage.rotate

c = np.zeros([10, 10]).astype(int)
for i in range(10):for j in range(10):c[i][j] = 10 * i + j
d = ndimage.rotate(c, 60, reshape=False)
print(c)
print(d)

上述代码就是将c这个二维数组（矩阵）逆时针旋转60度得到d，下面展示c和d的内容：

[[ 0  1  2  3  4  5  6  7  8  9][10 11 12 13 14 15 16 17 18 19][20 21 22 23 24 25 26 27 28 29][30 31 32 33 34 35 36 37 38 39][40 41 42 43 44 45 46 47 48 49][50 51 52 53 54 55 56 57 58 59][60 61 62 63 64 65 66 67 68 69][70 71 72 73 74 75 76 77 78 79][80 81 82 83 84 85 86 87 88 89][90 91 92 93 94 95 96 97 98 99]]
[[ 0  0  7 17 27 36  0  0  0  0][ 0  0 11 22 30 40 49 58  0  0][ 0  6 16 25 35 44 53 62 71 80][ 0 10 21 29 39 48 57 66 75 84][ 5 16 24 34 43 52 61 70 79 89][10 20 29 38 47 56 65 75 83 94][15 24 33 42 51 60 70 78 89  0][19 28 37 46 55 64 74 83 93  0][ 0  0 41 50 59 69 77 88  0  0][ 0  0  0  0 63 72 82 92  0  0]]

二、直接反投影

# 直接反投影
from scipy import ndimage
import numpy as np
import matplotlib.pyplot as pltIMG_SIZE = 512def IRadonTransform(sinogram, views):detectors = IMG_SIZEorigin = np.zeros((views, detectors, detectors))for i in range(views):projectionValue = sinogram[:, i]projectionValueExpandDim = np.expand_dims(projectionValue, axis=0)projectionValueRepeat = projectionValueExpandDim.repeat(IMG_SIZE, axis=0)origin[i] = ndimage.rotate(projectionValueRepeat, i*180/views, reshape=False).astype(np.float64)iradon = np.sum(origin, axis=0)return iradon# image size是512*512，探测器单元数量为512
# radon是256*512的sinogram，radon1是512*512的sinogram
# 即正投影时分别采用了256和512个views
img = IRadonTransform(radon, 256)
img1 = IRadonTransform(radon1, 512)plt.subplot(1, 2, 1)
plt.imshow(img, cmap='gray')
plt.subplot(1, 2, 2)
plt.imshow(img1, cmap='gray')
plt.show()

radon和radon1都是上一篇文章中的代码获取的sinogram，view分别是256和512，进行反投影的结果如下：

暂时还没有搞懂的是，为什么view的数量差那么多（一个是256一个是512），但是反投影的结果几乎一样。

三、先滤波再反投影

# 先滤波再反投影
import numpy as np
from scipy import ndimage
from scipy.signal import convolve
import matplotlib.pyplot as plt#两种滤波器的实现
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 filterSLIMG_SIZE = 512def IRadonTransform(sinogram, views):detectors = IMG_SIZEorigin = np.zeros((views, detectors, detectors))filter = SLFilter(detectors, 1)for i in range(views):projectionValue = sinogram[:, i]projectionValueFiltered = convolve(filter, projectionValue, "same")projectionValueExpandDim = np.expand_dims(projectionValueFiltered, axis=0)projectionValueRepeat = projectionValueExpandDim.repeat(IMG_SIZE, axis=0)origin[i] = ndimage.rotate(projectionValueRepeat, i*180/views, reshape=False).astype(np.float64)iradon = np.sum(origin, axis=0)return iradon# image size是512*512，探测器单元数量为512
# radon是256*512的sinogram，radon1是512*512的sinogram
# 即正投影时分别采用了256和512个views
FBP_img = IRadonTransform(radon, 256)
FBP_img1 = IRadonTransform(radon1, 512)plt.subplot(1, 2, 1)
plt.imshow(FBP_img, cmap='gray')
plt.subplot(1, 2, 2)
plt.imshow(FBP_img1, cmap='gray')
plt.show()

观察运行结果：

可以发现，view越多，条状伪影的数量越少。但是没搞懂的是，反投影既然是对不同view的投影数据的累加，那么按理来说view越多，应该像素数值越大。但实际上是view越多，反投影的结果越精确，这是为什么，资料有限，暂时没有搞清楚。

四、小结

这篇和上一篇文章是利用体模数据对CT image的正投影（获得sinogram）和反投影（通过sinogram获得CT image）进行了模拟，大概明白了CT成像重建的原理。但实际的原理肯定更为复杂，这就要靠之后继续学习了。

编程记录——研究一下python对shepp_logan体模数据实现iradon变换相关推荐

OpenCV函数应用：基于二值图像的三种孔洞填充方法记录（附python，C++代码）
系列文章目录函数系列: OpenCV函数简记_第一章数字图像的基本概念(邻域,连通,色彩空间) OpenCV函数简记_第二章数字图像的基本操作(图像读写,图像像素获取,图像ROI获取,图像混合,图形 ...
小白的第一本python书_读书笔记：编程小白的第一本python入门书
书名:编程小白的第一本python入门书作者:侯爵出版社/出处:图灵社区年份:2016年封面: 感想: 本书短小精悍,精华部分在于给编程小白打了鸡血的同时输出了一种"高效学习法的思想 ...
没学过编程可以自学python吗-完全没学过编程的人学习 Python前应该掌握些什么？...
在众多高大上的自学指导中,尝试做一股清流,把要讲清楚的都讲清楚,除了一堆资料之外,你能在学之前就有一个非常明显的结果倾向. 本文以<小白带你学Python>为内容方向,试图在繁杂的信息里, ...
python编程在哪里写-python入门该从哪里开始？
相信对于每个人而言,知道编程和学习编程这件事,出发点是不同的.汤哥在北京接触编程的时间是2013年,那个时候还在一个二线城市上大学,还没有这么多各种融资,各种互联网创业的氛围,大家想的更多的是一些线下 ...
python风变编程是骗局吗-风变编程：花时间学Python，是对自己未来最好的投资
风变编程:花时间学Python,是对自己未来最好的投资 2020年07月16日 16:00作者:网络编辑:王动分享谷歌研究主任Peter Norvig曾说:从一开始,Python就一直是谷歌的重要 ...
python编程300集免费-python 300本电子书合集
链接: https://pan.baidu.com/s/1CNlB35ASnDNlUGNCZJbiAA 提取码: fxig Q群:592857363 更多所在数据科学速查表零起点Python机器学 ...
python编程输入标准-揭秘python编程技巧
揭秘python编程技巧一.python的标准输入和输出[root@133 wc]# vim stdin.py #!/usr/bin/python #encoding:utf-8 import sy ...
从零开始学习python编程-如何从零开始学python？
在众多高大上的自学指导中,尝试做一股清流,把要讲清楚的都讲清楚,除了一堆资料之外,你能在学之前就有一个非常明显的结果倾向. 本文以<小白带你学Python>为内容方向,试图在繁杂的信息里, ...
python编程入门-Python编程入门经典pdf(Python编程入门教程) 高清中文版
Python编程入门经典pdf(Python编程入门教程)下载.Python编程入门经典pdf高清版帮助各位更好的进行Python编程的学习以及理解,最经典的课题,最深入的概念,让你在Python编程 ...
python趣味编程100例-儿童Python趣味编程课程
儿童Python趣味编程课程南京杜恩培训隶属于南京卡尔威特教育咨询有限公司,秉承"以人为本"的办学宗旨,致力于给每一位前来学习的学员专业优质的服务.中心成立于2002年,在南京已 ...

编程记录——研究一下python对shepp_logan体模数据实现iradon变换

目录

一、补充

二、直接反投影

三、先滤波再反投影

四、小结

编程记录——研究一下python对shepp_logan体模数据实现iradon变换相关推荐

最新文章

热门文章