Python实现ART重建算法_scipy中的art

作者：你好赵伟 | 2024-07-13 23:59:56

踩

scipy中的art

ART算法的简介

大家好，半个多月之前，我介绍了Radon变换和直接反投影以及滤波反投影的算法，现在向大家介绍一下ART算法，这是另一种CT图像重建的算法，同时给出Python实现。下面先简单地介绍一下投影矩阵的生成和ART算法的数学基础。
投影矩阵是代数重建算法的基础，它将投影数据和断层图像联系了起来，投影矩阵的计算方法也将影响重建图像的质量，投影矩阵的模型可以分为以下几种：

把射束看为是宽度为0，间距为 $\delta$ 的一系列直线，将尺寸为 $N * N$ 的断层的左上角的像素记为1号像素，右下角记为第 $N * N$ 号像素，如果第i条射线通过了j号像素，则记 $\omega_{ij}$ 为1，否则记为0，这样就形成了一个投影矩阵
$\begin{cases} 1, 第i条射束通过了第j号像素,\\ 0, 其他. \end{cases}$
把射束看作是宽度为零，间距为 $\delta$ 的一系列直线，将尺寸为 $N * N$ 的断层的左上角的像素记为1号像素，右下角记为第 $N * N$ 号像素，如果第i条射线通过了j号像素，则记 $\omega_{ij}$ 为1，否则记为0，同时考虑射束在它经过的像素中走过的距离。
把射束看作是一个宽度为 $l$ 的一条带，将尺寸为 $N * N$ 的断层的左上角的像素记为1号像素，右下角记为第 $N * N$ 号格点，如果第i条射束的带通过了第j号像素的中心点，则记 $\omega_{ij}$ 为1，否则为0
$\begin{cases} 1, 第i条射束通过了第j号像素的中心点,\\ 0, 其他. \end{cases}$
把射束看作是一个宽度为 $l$ 的一条带，将尺寸为 $N * N$ 的断层的左上角的像素记为1号像素，右下角记为第 $N * N$ 号格点，计算射束带与它所经过的像素点相交的面积与像素的面积之比，将比值记为 $\omega_{ij}$ ，形成投影矩阵。

其中模型1、模型3较为简单，但是容易受噪声影响，而模型4 计算量较大，所以后面的源码中实现的是模型2。

Algebraic Reconstruction Technique (ART) 代数重建算法具有可以利用不完全的信息进行重建（不需要将180 $\degree$ 分为几百上千个步来扫描，而FBP算法需要的投影数据较多），重建结果不易受噪声影响等优点，但是重建耗费的时间长。这个算法其实就是解一个方程，射束组成一个一维数组，记为 $p_{1}, p_{2}, ..., p_{m}]$ ，其中 $m$ 为穿过物体的x射线的条数，它等于投影角度与每个投影角度上射束的条目的乘积，则有一个方程：

$\begin{cases} \omega_{11} f_1 + \omega_{12}f_2+...+\omega_{1N} f_N = p_{1}\\ \omega_{21} f_1 + \omega_{22}f_2+...+\omega_{2N} f_N = p_{2}\\ ...\\ \omega_{m1} f_1 + \omega_{m2}f_2+...+\omega_{mN} f_N = p_{m} \end{cases}$
一共N个未知数（ $f_{1}, f_{2}, ..., f_{N}$ ），M个方程，其中的 $\omega_{ij}$ 就是前面介绍的投影矩阵的元素。我们只要进行以下形式的迭代：
$f_{j}^{k+1} = f_{j}^{k} + \lambda \frac{p_{i}-\sum_{n=1}^{N} \omega_{in}f_{n}^{k}}{\sum_{n=1}^{N}\omega_{in}^{2}}\omega_{ij}, \quad j=1, 2, ..., N$
进行一次迭代就是一次投影和反投影操作，我们可以进行多次迭代，直到结果收敛。

Python源码

#作者：CSDN用户“宋体的微软雅黑（hsyxxyg）”
#时间：2020年6月18日
#脚本任务：生成投影系统矩阵，并利用此矩阵进行ART重建。

import numpy as np
import pandas as pd
from scipy import ndimage
import imageio
import matplotlib.pyplot as plt


def CalSystemMatrix(theta, pictureSize, projectionNum, delta):
    squareChannels = np.power(pictureSize, 2)
    totalPorjectionNum = len(theta) * projectionNum
    gridNum = np.zeros((totalPorjectionNum, 2 * pictureSize))
    gridLen = np.zeros((totalPorjectionNum, 2 * pictureSize))
    t = np.arange(-(pictureSize - 1) / 2, (pictureSize - 1) / 2+1)
    for loop1 in range(len(theta)):
        for loop2 in range(pictureSize):
            u = np.zeros((2 * pictureSize))
            v = np.zeros((2 * pictureSize))
            th = theta[loop1]
            if th == 90:
            #如果计算的结果超出了网格的范围，则立刻开始计算下一个射束
                if ((t[loop2] >= pictureSize / 2 * delta) or (t[loop2] <= - pictureSize / 2 * delta)):
                    continue
                kout = pictureSize * np.ceil(pictureSize/2 - t[loop2]/delta)
                kk = np.arange(kout - (pictureSize -1 ), kout+1)
                u[0:pictureSize] = kk
                v[0:pictureSize] = np.ones(pictureSize) * delta
                
            elif th==0:
                if (t[loop2] >= pictureSize / 2 * delta) or (t[loop2] <= -pictureSize / 2 * delta):
                    continue
                kin = np.ceil(pictureSize/2 + t[loop2] / delta)
                kk = np.arange(kin, (kin + pictureSize * pictureSize), step=pictureSize)
                u[0:pictureSize] = kk
                v[0:pictureSize] = np.ones(pictureSize) * delta

            else:
                if th>90:
                    th_temp = th - 90
                elif th<90:
                    th_temp = 90 - th
                th_temp = th_temp * np.pi / 180
                #计算束线的斜率和截距
                b = t / np.cos(th_temp)
                m = np.tan(th_temp)
                y1d = -(pictureSize / 2) * delta * m + b[loop2]
                y2d = (pictureSize / 2) * delta * m + b[loop2]
                #if (y1d < -pictureSize / 2 * delta and y2d < -pictureSize/2 * delta) or (y1d > pictureSize / 2 * delta and y2d > -pictureSize / 2 * delta):
                if (y1d < -pictureSize / 2 * delta and y2d < -pictureSize/2 * delta) or (y1d > pictureSize / 2 * delta and y2d > pictureSize / 2 * delta):
                    continue
                if (y1d <= pictureSize / 2 * delta and  y1d >= -pictureSize/2 * delta and y2d > pictureSize / 2 * delta):
                    yin = y1d
                    d1 = yin - np.floor(yin / delta) * delta
                    kin = pictureSize * np.floor(pictureSize / 2 - yin / delta) + 1
                    yout = pictureSize / 2 * delta
                    xout = (yout - b[loop2]) / m
                    kout = np.ceil(xout / delta) + pictureSize / 2

                elif (y1d <= pictureSize/2 * delta and y1d >= -pictureSize/2 * delta and y2d >= -pictureSize/2 * delta and y2d < pictureSize/2 * delta):
                    yin = y1d
                    d1 = yin - np.floor(yin/delta) * delta
                    kin = pictureSize * np.floor(pictureSize / 2 - yin / delta) + 1
                    yout = y2d
                    #1:
                    #2:xout = (yout - b[loop2]) / m
                    kout = pictureSize * np.floor(pictureSize/2 - yout/delta) + pictureSize

                elif (y1d < - pictureSize / 2 * delta and y2d > pictureSize / 2 * delta):
                    yin = - pictureSize / 2  * delta
                    xin = (yin - b[loop2]) / m
                    d1 = pictureSize / 2 * delta + np.floor(xin / delta) * delta * m + b[loop2]
                    kin = pictureSize * (pictureSize - 1) + pictureSize / 2 + np.ceil(xin / delta)
                    yout = pictureSize / 2 * delta
                    #error: xout = (yout / b[loop2])/m
                    xout = (yout - b[loop2]) / m
                    kout = np.ceil(xout / delta) + pictureSize / 2

                elif (y1d < - pictureSize / 2 * delta and y2d >= -pictureSize / 2 * delta and y2d < pictureSize / 2 * delta):
                    yin = -pictureSize / 2 * delta 
                    xin = (yin - b[loop2]) / m
                    d1 = pictureSize / 2 * delta + np.floor(xin / delta) * delta * m + b[loop2]
                    kin = pictureSize * (pictureSize - 1) + pictureSize / 2 + np.ceil(xin / delta)
                    yout = y2d
                    kout = pictureSize * np.floor(pictureSize / 2 - yout / delta) + pictureSize
                else:
                    continue
                #计算第i条射束穿过的像素的编号和长度
                k = kin
                c = 0
                d2 = d1 + m * delta
                while k >= 1 and k <= squareChannels:
                    if d1 >= 0 and d2 > delta:
                        u[c] = k
                        v[c] = (delta - d1) * np.sqrt(np.power(m, 2) + 1) / m
                        if k > pictureSize and k != kout:
                            k = k - pictureSize
                            d1 = d1 - delta
                            d2 = d1 + m * delta
                        else:
                            break
                    elif d1 >= 0 and d2 == delta:
                        u[c] = k
                        v[c] = delta * np.sqrt(np.power(m, 2) + 1)
                        if k>pictureSize and k != kout:
                            k = k - pictureSize + 1
                            d1 = 0
                            d2 = d1 + m * delta
                        else:
                            break
                    elif d1 >= 0 and d2 < delta:
                        u[c] = k
                        v[c] = delta * np.sqrt(np.power(m, 2) + 1)
                        if k!=kout:
                            k = k + 1
                            d1 = d2
                            d2 = d1 + m * delta
                        else:
                            break
                    elif d1 <= 0 and d2 >= 0 and d2 <= delta:
                        u[c] = k
                        v[c] = d2 * np.sqrt(np.power(m, 2) + 1) / m
                        if k != kout:
                            k = k + 1
                            d1 = d2
                            d2 = d1 + m * delta
                        else:
                            break
                    elif d1 <= 0 and d2 > delta:
                        u[c] = k
                        v[c] = delta * np.sqrt(np.power(m, 2) + 1) / m
                        if k > pictureSize and k != kout:
                            k = k - pictureSize
                            #k = k + 1
                            d1 = -delta + d1
                            d2 = d1 + m * delta
                        else:
                            break
                    else:
                        print(d1, d2, "error!!!")
                    c = c + 1
                #如果投影角度小于90度，应该利用投影射线关于y轴的对称性计算出权重因子向量。
                if th < 90:
                    u_temp = np.zeros(2 * pictureSize)
                    if u.any() == 0:
                        continue
                    indexULTZero = np.where(u>0)
                    for innerloop in range(len(u[indexULTZero])):
                        r = np.mod(u[innerloop], pictureSize)
                        if r == 0:
                            u_temp[innerloop] = u[innerloop] - pictureSize
                        else:
                            u_temp[innerloop] = u[innerloop] - 2 * r + pictureSize
                    u = u_temp
            gridNum[loop1 * projectionNum + loop2, :] = u
            gridLen[loop1 * projectionNum + loop2, :] = v
    return gridNum, gridLen

def DiscreteRadonTransform(image, theta):
    projectionNum = len(image[0])
    thetaLen = len(theta)
    radontansformRes = np.zeros((projectionNum, thetaLen), dtype='float64')
    for s in range(len(theta)):
        rotation = ndimage.rotate(image, -theta[s], reshape=False).astype('float64')
        radontansformRes[:, s] = sum(rotation)
    return radontansformRes
#适用于灰度图
image = imageio.imread("shepplogan.jpg").astype(np.float64)
#如果你的图片的色彩模式是RGB或者RGBA请使用以下语句
#如果使用上述imageio.imread，会出现数组维度不匹配的错误。
#image = cv2.imread("shepplogan.jpg", cv2.IMREAD_GRAYSCALE)

theta = np.linspace(0, 180, 60, dtype=np.float64)
#使用离散Radon变换获取投影值
projectionValue = DiscreteRadonTransform(image, theta)

#定义用到一些参数：旋转角度的矩阵，探测器的道数，图片尺寸，平移步长，最大迭代次数，驰豫因子
projectionNum = np.int64(256)
pictureSize = np.int64(256)
pictureSizeSquare = pictureSize * pictureSize
delta = np.int64(1)
irt_Num = np.int64(20)
lam = np.float64(0.25)

#计算投影矩阵
gridNum, gridLen = CalSystemMatrix(theta, pictureSize, projectionNum, delta)

dfgridNum = pd.DataFrame(gridNum)
dfgridLen = pd.DataFrame(gridLen)
#可以将系统矩阵存储到文件中,以后直接使用。
dfgridNum.to_csv("gridNum.csv", header=False, index=False)
dfgridLen.to_csv("gridLen.csv", header=False, index=False)

#gridNum = np.array(pd.read_csv("gridNum1.csv"), header=None)
#gridLen = np.array(pd.read_csv("gridLen1.csv"), header=None)
#存储重建获得的图像的矩阵
F = np.zeros((pictureSize*pictureSize, ))

irt_Num = 10
lam = 0.25
c = 0

#开始迭代过程
while(c < irt_Num):
    for loop1 in range(len(theta)):
        for loop2 in range(pictureSize):
            u = gridNum[loop1 * pictureSize + loop2, :]
            v = gridLen[loop1 * pictureSize + loop2, :]
            if u.any() == 0:
                continue
            w = np.zeros(pictureSizeSquare, dtype=np.float64)
            uLargerThanZero = np.where(u > 0)
            w[u[uLargerThanZero].astype(np.int64)-1] = v[uLargerThanZero]
            PP = w.dot(F)
            C = (projectionValue[loop2, loop1] - PP) / sum(np.power(w, 2)) * w.conj()
            F = F + lam * C
    F[np.where(F < 0)] = 0
    c = c + 1

F = F.reshape(pictureSize, pictureSize).conj()

#绘制Sinogram图和重建的结果
plt.subplot(1, 2, 1)
plt.imshow(projectionValue, cmap="gray")

plt.subplot(1, 2, 2)
plt.imshow(F, cmap="gray")
plt.savefig("modify90Sample.png", 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
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
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231

代码运行的结果：
上述代码运行的结果

上述代码运行得到的结果，左侧是投影值，右侧是重建的结果。

写在最后：上面的代码已经是可以运行的程序了，但是计算时间有点长，在我的CPU为Ryzen 4800u的笔记本上需要27分钟左右的时间（包括生成投影矩阵和迭代过程）。当然，我们可以使用numba对迭代过程中的循环进行加速，我大致看了一下，加速之后代码的运行时间缩短到了4分钟左右。大家可以自行尝试一下使用numba对Python的循环过程进行提速，还挺简单的。

本文内容由网友自发贡献，转载请注明出处：https://www.wpsshop.cn/w/你好赵伟/article/detail/822218