赞
踩
所有图像处理重建可在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()

Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。