Python实现ART重建算法

宋体的微软雅黑

已于 2024-01-06 13:52:53 修改

阅读量7.7k

点赞数 44

文章标签： python

于 2020-06-26 21:28:30 首次发布

本文链接：https://blog.csdn.net/hsyxxyg/article/details/106962790

版权

ART算法的简介

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

把射束看为是宽度为0，间距为 $\delta$ 的一系列直线，将尺寸为 $N * N$ 的断层的左上角的像素记为1号像素，右下角记为第 $N * N$ 号像素，如果第i条射线通过了j号像素，则记 $\omega_{ij}$ 为1，否则记为0，这样就形成了一个投影矩阵
$\omega_{ij} = \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
$\omega_{ij} = \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()

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

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

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