1.ART方法：

1.代码实现：

import numpy as np
import matplotlib.pyplot as plt
import cv2
from scipy import ndimage
import timedef projection(image, theta):"""计算投影值:param image: 原始图像:param theta: 射束旋转角度:return: 投影值矩阵projectionNum: 射束个数thetaNum: 角度个数randontansform: 拉登变换结果"""print('正在计算投影...')projectionNum = len(image[0])thetaNum = len(theta)radontansform = np.zeros((projectionNum, thetaNum), dtype='float64')for i in range(len(theta)):# 进行离散拉登变换rotation = ndimage.rotate(image, -theta[i], reshape=False).astype('float64')radontansform[:, i] = sum(rotation)return radontansformdef systeMartrix(theta, size, num, delta):"""计算系统矩阵:param theta: 射束旋转角度:param size: 图片尺寸:param num: 射束条数:param delat: 网格边长:return: gridNum：穿过网格编号，gridLen：穿过网格长度"""print('正在计算系统矩阵...')start = time.perf_counter()functionNum = len(theta) * num  # 方程个数，即为系统矩阵的行数# 射线最多穿过2*size个方格gridNum = np.zeros((functionNum, 2 * size))  # 系统矩阵：编号gridLen = np.zeros((functionNum, 2 * size))  # 系统矩阵：长度N = np.arange(-(size - 1) / 2, (size - 1) / 2 + 1)  # 射束for loop1 in range(len(theta)):th = theta[loop1]  # 射束角度for loop2 in range(size):u = np.zeros((2 * size))  # 编号v = np.zeros((2 * size))  # 长度# 垂直入射if th == 0:# 射束未穿过图像时if (N[loop2] >= size / 2 * delta) or (N[loop2] <= -size / 2 * delta):continue# 入射网格编号kin = np.ceil(size / 2 + N[loop2] / delta)# 穿过网格编号kk = np.arange(kin, (kin + size * size), step=size)u[0:size] = kkv[0:size] = np.ones(size) * delta# 平行入射elif th == 90:if (N[loop2] >= size / 2 * delta) or (N[loop2] <= -size / 2 * delta):continue# 出射网格编号kout = size * np.ceil(size / 2 - N[loop2] / delta)kk = np.arange(kout - size + 1, kout + 1)u[0:size] = kkv[0:size] = np.ones(size) * deltaelse:# phi为射束与x轴所夹锐角if th > 90:phi = th - 90elif th < 90:phi = 90 - th# 角度值换算为弧度制phi = phi * np.pi / 180# 截距b = N / np.cos(phi)# 斜率m = np.tan(phi)# 入射点纵坐标y1 = -(size / 2) * delta * m + b[loop2]# 出射点纵坐标y2 = (size / 2) * delta * m + b[loop2]# 射束未穿过图像if (y1 < -size / 2 * delta and y2 < -size / 2 * delta) or (y1 > size / 2 * delta and y2 > size / 2 * delta):continue# 穿过a、b边（左侧和上侧）if (y1 <= size / 2 * delta and y1 >= -size / 2 * delta and y2 > size / 2 * delta):"""(xin,yin): 入射点坐标(xout,yout): 出射点坐标kin,kout: 入射格子标号，出射格子编号d1: 入射格子左下角与入射射束距离"""yin = y1yout = size / 2 * delta# xin = -size / 2 * deltaxout = (yout - b[loop2]) / mkin = size * np.floor(size / 2 - yin / delta) + 1kout = np.ceil(xout / delta) + size / 2d1 = yin - np.floor(yin / delta) * delta# 穿过a、c边（左侧和右侧）elif (y1 <= size / 2 * delta and y1 >= -size / 2 * delta and y2 >= -size / 2 * delta and y2 < size / 2 * delta):# xin = -size / 2 * delta# xout = size / 2 * deltayin = y1yout = y2kin = size * np.floor(size / 2 - yin / delta) + 1kout = size * np.floor(size / 2 - yout / delta) + sized1 = yin - np.floor(yin / delta) * delta# 穿过d、b边（下侧和上侧）elif (y1 < - size / 2 * delta and y2 > size / 2 * delta):yin = - size / 2 * deltayout = size / 2 * deltaxin = (yin - b[loop2]) / mxout = (yout - b[loop2]) / mkin = size * (size - 1) + size / 2 + np.ceil(xin / delta)kout = np.ceil(xout / delta) + size / 2d1 = size / 2 * delta + np.floor(xin / delta) * delta * m + b[loop2]# 穿过d、c边（下侧和右侧）elif (y1 < - size / 2 * delta and y2 >= -size / 2 * delta and y2 < size / 2 * delta):yin = -size / 2 * deltayout = y2xin = (yin - b[loop2]) / m# xout = size / 2 * deltakin = size * (size - 1) + size / 2 + np.ceil(xin / delta)kout = size * np.floor(size / 2 - yout / delta) + sized1 = size / 2 * delta + np.floor(xin / delta) * delta * m + b[loop2]else:continue# 计算穿过的格子编号和长度"""k: 射线穿过的格子编号c: 穿过格子的序号d2: 穿过的格子的右侧与该方格右下角顶点的距离"""k = kin  # 入射的格子即为穿过的第一个格子c = 0  # c为穿过格子的序号d2 = d1 + m * delta  # 与方格右侧交点# 当方格数在1~n^2内迭代计算while k >= 1 and k <= np.power(size, 2):"""根据射线与方格的左右两侧的交点关系，来确定穿过方格的六种情况。在每种情况中，存入穿过的方格编号，穿过方格的射线长度。若该方格是最后一个穿过的方格，则停止迭代；若不是最后一个方格，则计算下一个穿过的方格的编号、左右边与射线的交点。"""if d1 >= 0 and d2 > delta:u[c] = k  # 穿过方格的编号v[c] = (delta - d1) * np.sqrt(np.power(m, 2) + 1) / m  # 穿过方格的射线长度if k > size and k != kout:  # 若不是最后一个方格k -= size  # 下一个方格编号d1 -= delta  # 下一个方格左侧交点d2 = delta * m + d1  # 下一个方格右侧交点else:  # 若是最后一个方格则直接跳出循环breakelif d1 >= 0 and d2 == delta:u[c] = kv[c] = delta * np.sqrt(np.power(m, 2) + 1)if k > size and k != kout:k -= size + 1d1 = 0d2 = d1 + m * deltaelse:breakelif d1 >= 0 and d2 < delta:u[c] = kv[c] = delta * np.sqrt(np.power(m, 2) + 1)if k != kout:k += 1d1 = d2d2 = d1 + m * deltaelse:breakelif d1 <= 0 and d2 >= 0 and d2 <= delta:u[c] = kv[c] = d2 * np.sqrt(np.power(m, 2) + 1) / mif k != kout:k += 1d1 = d2d2 = d1 + m * deltaelse:breakelif d1 <= 0 and d2 > delta:u[c] = kv[c] = delta * np.sqrt(np.power(m, 2) + 1) / mif k > size and k != kout:k -= sized1 -= deltad2 = d1 + m * deltaelse:breakelif d1 <= 0 and d2 == delta:u[c] = kv[c] = d2 * np.sqrt(np.power(m, 2) + 1) / mif k > size and k != kout:k -= size + 1d1 = 0d2 = delta * melse:breakelse:print(d1, d2, "数据错误！")c += 1# 当射线斜率为负数时，利用对称性进行计算if th < 90:u_temp = np.zeros(2 * size)# 排除掉未穿过图像的射束if u.any() == 0:continue# 射束穿过的格子的编号indexMTZero = np.where(u > 0)# 利用对称性求得穿过网格编号for loop in range(len(u[indexMTZero])):"""计算穿过的方格编号，利用方格编号与边长的取余的关系，得到对称的射线穿过的方格编号。"""r = np.mod(u[loop], size)if r == 0:u_temp[loop] = u[loop] - size + 1else:u_temp[loop] = u[loop] - 2 * r + size + 1u = u_temp# 方格编号gridNum[loop1 * num + loop2, :] = u# 穿过方格的射线长度gridLen[loop1 * num + loop2, :] = vend = time.perf_counter()print('系统矩阵计算耗时：%s 秒。'%(end - start))return gridNum, gridLendef iteration(theta, size, gridNum, gridLen, F, ite_num):"""按照公式迭代重建:param theta: 旋转角度:param size: 图像边长:param gridNum: 射线穿过方格编号:param gridLen: 射线穿过方格长度:param F: 重建后图像:return: 重建后图像"""print('正在进行迭代...')start = time.perf_counter()c = 0  # 迭代计数while (c < ite_num):for loop1 in range(len(theta)):  # 在角度theta下for loop2 in range(size):  # 第loop2条射线u = gridNum[loop1 * size + loop2, :]v = gridLen[loop1 * size + loop2, :]if u.any() == 0:  # 若射线未穿过图像，则直接计算下一条射线continue# 本条射线对应的行向量w = np.zeros(sizeSquare, dtype=np.float64)# 本条射线穿过的网格编号uLargerThanZero = np.where(u > 0)# 本条射线穿过的网格长度w[u[uLargerThanZero].astype(np.int64) - 1] = v[uLargerThanZero]# 计算估计投影值PP = w.dot(F)# 计算实际投影与估计投影的误差error = projectionValue[loop2, loop1] - PP# 求修正值C = error / sum(np.power(w, 2)) * w.conj()# 进行修正F = F + lam * CF[np.where(F < 0)] = 0c = c + 1F = F.reshape(size, size).conj()end = time.perf_counter()print('迭代耗时：%s 秒。'%(end - start))return Fif __name__ == '__main__':image = cv2.imread("shepp-logan(256).png", cv2.IMREAD_GRAYSCALE)"""theta: 射束旋转角度num: 射束条数size: 图片尺寸delta: 网格边长ite_num: 迭代次数lam: 松弛因子"""theta = np.linspace(0, 180, 60, dtype=np.float64)num = np.int64(256)size = np.int64(256)delta = np.int64(1)ite_num = np.int64(10)lam = np.float64(.25)sizeSquare = size * size# 计算投影值projectionValue = projection(image, theta)# 计算系统矩阵gridNum, gridLen = systeMartrix(theta, size, num, delta)# 重建后图像矩阵F = np.zeros((size * size,))# 迭代法重建reconImage = iteration(theta, size, gridNum, gridLen, F, ite_num)print('共迭代%d次。'%ite_num)# 绘制原始图像、重建图像、误差图像plt.subplot(1, 3, 1)plt.imshow(image, cmap='gray')plt.title('Original Image')plt.subplot(1, 3, 2)plt.imshow(reconImage, cmap='gray')plt.title('Reconstruction Image')plt.subplot(1, 3, 3)plt.imshow(reconImage - image, cmap='gray')plt.title('Error Image')plt.savefig("ART.png", cmap='gray')plt.show()

2.耗时：

3.重建结果：

2.SART方法：

1.代码实现：

对迭代过程进行修改：

def iteration(theta, size, gridNum, gridLen, F, ite_num, num, sizeSquare, projectionValue, lam):"""按照公式迭代重建:param theta: 旋转角度:param size: 图像边长:param gridNum: 射线穿过方格编号:param gridLen: 射线穿过方格长度:param F: 重建后图像:param sizeSquare: 图像边长平方:param projectionValue: 投影值:param lam: 松弛因子:return: 重建后图像"""print('正在进行迭代...')start = time.perf_counter()c = 0  # 迭代计数while (c < ite_num):for loop1 in range(len(theta)):  # 在角度theta下u = gridNum[loop1 * num : (loop1 + 1) * num, :]  # 射线穿过网格编号v = gridLen[loop1 * num : (loop1 + 1) * num, :]  # 射线穿过网格长度w = np.zeros((num, sizeSquare))for loop2 in range(num):  # 对第loop2条射线if u[loop2, :].any() == 0:  # 若射线未穿过图像，则直接计算下一条射线continuefor loop3 in range(2 * size):m = u[loop2, loop3]if ((m > 0) and (m <= sizeSquare)):w[loop2, np.int64(m)-1] = v[loop2, loop3]sumcol = np.sum(w, axis=0).reshape(-1, 1)sumrow = np.sum(w, axis=1).reshape(-1, 1)ind1 = np.where(sumrow > 0)corr = np.zeros((num, 1))err = projectionValue[:, loop1].reshape((256, 1)) - w.dot(F)corr[ind1] = err[ind1] * (1 / sumrow[ind1])backPro = (w.T).dot(corr)ind2 = np.where(sumcol > 0)delta = np.zeros((sizeSquare, 1))delta[ind2] = backPro[ind2] / sumcol[ind2]F += lam * deltaF[np.where(F < 0)] = 0c += 1F = F.reshape(size, size)end = time.perf_counter()print('迭代耗时：%s 秒。' % (end - start))return F

注意此时在main中创建F时应改为：

    # 重建后图像矩阵F = np.zeros((sizeSquare, 1))

2.耗时：

3.重建结果：

ART与SART代数重建算法相关推荐

代数重建算法（ART）详解
迭代图像重建的方法可分为代数迭代法和统计迭代法,代数迭代法以代数重建算法(Algebraic Reconstruction Technique,ART)为代表.ART适合于不完全投影数据的图像重建,其 ...
迭代重建算法中投影矩阵的计算
在前面学习的重建算法都是属于解析法,它是以Radon变换为理论基础,首先对投影数据在连续域进行一些处理,然后将其离散化进行重建.接下来学习到的是图像重建的迭代算法,该算法的主要思路是:从一幅假设的初始 ...
口腔ct重建服务器原理,口腔CT成像的迭代重建算法研究
摘要: 自CT技术诞生以来,它就被认为是二十世纪影响人类发展的十大技术之一.其应用更是涉及了多个领域,包括医学,生物,工业,安全等等.特别是在医学方面,CT更是成为了医学检测的主要手段.而口腔CT(D ...
matlab_医学CT重建 ART，SART算法
一,代码下载 matlab_医学CT重建ART,SART算法-自然语言处理文档类资源-CSDN下载二.ART算法 1.基本思想 ART迭代算法的基本思想是先将连续的图像离散化,再采用CT成像的离散模 ...
二维温度场matlab编程,二维温度场重建算法（价钱可议）
即病态方程的求解,具体问题如下: (1)物理模型矩形区域的物理模型,即辐射传热定律: 可转化为:E=UT 其中:E表示n个壁面单元辐射接收装置各自所接收的辐射能组成的集合. U被称为敏感矩阵,其与第 ...
研究 | CT图像迭代重建算法研究进展
上次讲到我实现了一下代数迭代重建(ART),到周六开会的时候才大概了解了我的研究方向应该是统计迭代重建,一下子就把我给搞懵了.按照书上的说法,统计迭代法是在发射型CT(SPECT和PET)中应用广泛, ...
matlab重建算法stomp,压缩感知图像重建算法的研究现状及其展望.pdf
压缩感知图像重建算法的研究现状及其展望.pdf D陆33EE iE γ,理想想 [本文献信息}李然,干宗良,崔子冠,等.压缩感知图像重建算法的研究现状及其展望[J]. 电视技术,20日,37 (19) ...
汇总|3D人脸重建算法
作者:Tom Hardy Date:2019-12-30 来源:汇总|3D人脸重建算法
K单体型重建算法的研究
K单体型重建算法的研究王兆灿 [摘要]:随着新一代基因测序技术的飞速发展,以及单体型数据在人类遗传学等领域研究和应用的不断深入,对单体型数据的研究开始转向其他生物物种.由于测序技术的限制,通过生 ...

ART与SART代数重建算法

1.ART方法：

1.代码实现：

2.耗时：

3.重建结果：

2.SART方法：

1.代码实现：

2.耗时：

3.重建结果：

ART与SART代数重建算法相关推荐

最新文章

热门文章