python投影数据预处理与重建_投影重建图像python

所有图像处理重建可在https://github.com/Ruitao-chen/Medical-Computed-Tomography-CT-Imaging-Experiment找到，记得给我一个star星星

已有360个角度扫描的样本数据和一个角度的空场数据，且探测器的暗场影响忽略不计。对数据进行预处理。依次进行平场校正、坏点补偿和对数操作，最后利用滤波反投影重建算法进行断层重建。

平场校正：

将512*1的列向量求其元素的均值mean，再将列向量扩充为512*360的矩阵得到bright_mean

坏点补偿：

查看矩阵发现 100和220行数据异常，因此选择上下两侧的数据加权求均值进行修正。

I99=(I98*0.3+I99*0.7+I101*0.7+I102*0.3)*0.5

I220=(I218*0.3+I219*0.7+I221*0.7+I222 *0.3 )*0.5

对数运算：

I0 = max(max(I));

I=log(I/I0);

·  数据处理流程：

·  打开并读取投影数据（ParallelBeam_ProjectionDara.raw）和平场数据（ParallelBeam_DetectorFlatData.raw）。

进行平场校正：将投影数据除以平场数据，再乘以平场数据的均值。

进行坏点补偿：通过加权平均的方式，修复探测器上的坏点。

进行对数运算操作：对平场校正后的数据进行取对数操作。

使用滤波反投影重建算法：对处理后的数据进行反投影重建，得到最终的断层图像。

·  展示处理步骤的结果：

·  平场校正前后的投影数据对比。

坏点补偿前后的投影数据对比。

对数运算前后的投影数据对比。

最终的断层重建图像。

·  分析各个步骤对重建结果的影响：

·  平场校正：平场校正可以消除探测器的非均匀性，使投影数据更加准确。

坏点补偿：坏点补偿可以修复探测器上的异常值，提高数据的质量。

对数运算：对数运算可以压缩动态范围，增强低灰度值的细节信息。

滤波反投影重建算法：通过反投影重建算法，可以将投影数据逆向投影到空间中，得到样品的断层图像。

实验结果：

代码

import numpy as np
import matplotlib.pyplot as plt
 
from skimage.transform import iradon
from matplotlib.font_manager import FontProperties
fname = '/home/pai4090/anaconda3/lib/python3.11/site-packages/matplotlib/mpl-data/fonts/ttf/SimHei.ttf'
myfont = FontProperties(fname=fname)
 
 
 
filename = '/home/pai4090/断层成像/5/ParallelBeam_DetectorFlatData.raw'
fid = open(filename, 'rb')
A = np.fromfile(fid, dtype=np.ushort).reshape(1, 512)
 
filename = '/home/pai4090/断层成像/5/ParallelBeam_ProjectionDara.raw'
fid = open(filename, 'rb')
B = np.fromfile(fid, dtype=np.ushort).reshape(360,512)
 
A = A.T
B = B.T
 
plt.figure(1)
plt.imshow(B, cmap='gray')
plt.title('data')
plt.show()  
# plt.figure()
# plt.imshow(A, cmap='gray')
# plt.title('data')
# plt.show()  
 
I_bright_mean = np.mean(A)
 
# 平场校正
I_sample_dark = B
I_bright_dark = np.tile(A, (1, 360))
 
I = (I_sample_dark / I_bright_dark) * I_bright_mean
I_noflat = I_sample_dark.copy()
 
# 坏点补偿
I_nobad = I.copy()
I[99, :] = (I[97, :] * 0.3 + I[98, :] * 0.7 + I[100, :] * 0.7 + I[101, :] * 0.3) * 0.5
I[219, :] = (I[217, :] * 0.3 + I[218, :] * 0.7 + I[220, :] * 0.7 + I[221, :] * 0.3) * 0.5
 
I_noflat[99, :] = (I_noflat[97, :] * 0.3 + I_noflat[98, :] * 0.7 + I_noflat[100, :] * 0.7 + I_noflat[101, :] * 0.3) * 0.5
I_noflat[219, :] = (I_noflat[217, :] * 0.3 + I_noflat[218, :] * 0.7 + I_noflat[220, :] * 0.7 + I_noflat[221, :] * 0.3) * 0.5
 
# 对数操作
I_nolog = I.copy()
I0 = np.max(np.max(I))
I = np.log(I / I0)
 
I0 = np.max(np.max(I_noflat))
I_noflat = np.log(I_noflat / I0)
 
I0 = np.max(np.max(I_nobad))
I_nobad = np.log(I_nobad / I0)
 
# 反投影重建
data = iradon(I, theta=list(range(360)), filter_name='ramp', interpolation='linear')
noflat = iradon(I_noflat, theta=list(range(360)), filter_name='ramp', interpolation='linear')
nobad = iradon(I_nobad, theta=list(range(360)), filter_name='ramp', interpolation='linear')
nolog = iradon(I_nolog, theta=list(range(360)), filter_name='ramp', interpolation='linear')
 
#显示结果
plt.figure(2)
plt.imshow(data, cmap='gray')
plt.title('重建结果',fontproperties=myfont)
plt.figure(3)
plt.imshow(noflat, cmap='gray')
plt.title('无平场校正',fontproperties=myfont)
plt.figure(4)
plt.imshow(nobad, cmap='gray')
plt.title('无坏点补偿',fontproperties=myfont)
plt.figure(5)
plt.imshow(nolog, cmap='gray')
plt.title('无对数操作',fontproperties=myfont)
plt.show()