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

相关文章：编程记录——研究一下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)
1
2
3
4
5
6
7

上述代码就是将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]]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

二、直接反投影

# 直接反投影
from scipy import ndimage
import numpy as np
import matplotlib.pyplot as plt

IMG_SIZE = 512

def IRadonTransform(sinogram, views):
    detectors = IMG_SIZE
    origin = 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()
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

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] = 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

IMG_SIZE = 512

def IRadonTransform(sinogram, views):
    detectors = IMG_SIZE
    origin = 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()
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

观察运行结果：

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

四、小结

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