CT重建算法，迭代重建算法，ART算法的实现（亲测可运行）

#作者：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 + l