当前位置:   article > 正文

python投影数据预处理与重建_投影重建图像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)。

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

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

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

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

·  展示处理步骤的结果:

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

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

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

最终的断层重建图像。

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

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

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

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

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

实验结果:

代码

  1. import numpy as np
  2. import matplotlib.pyplot as plt
  3. from skimage.transform import iradon
  4. from matplotlib.font_manager import FontProperties
  5. fname = '/home/pai4090/anaconda3/lib/python3.11/site-packages/matplotlib/mpl-data/fonts/ttf/SimHei.ttf'
  6. myfont = FontProperties(fname=fname)
  7. filename = '/home/pai4090/断层成像/5/ParallelBeam_DetectorFlatData.raw'
  8. fid = open(filename, 'rb')
  9. A = np.fromfile(fid, dtype=np.ushort).reshape(1, 512)
  10. filename = '/home/pai4090/断层成像/5/ParallelBeam_ProjectionDara.raw'
  11. fid = open(filename, 'rb')
  12. B = np.fromfile(fid, dtype=np.ushort).reshape(360,512)
  13. A = A.T
  14. B = B.T
  15. plt.figure(1)
  16. plt.imshow(B, cmap='gray')
  17. plt.title('data')
  18. plt.show()
  19. # plt.figure()
  20. # plt.imshow(A, cmap='gray')
  21. # plt.title('data')
  22. # plt.show()
  23. I_bright_mean = np.mean(A)
  24. # 平场校正
  25. I_sample_dark = B
  26. I_bright_dark = np.tile(A, (1, 360))
  27. I = (I_sample_dark / I_bright_dark) * I_bright_mean
  28. I_noflat = I_sample_dark.copy()
  29. # 坏点补偿
  30. I_nobad = I.copy()
  31. I[99, :] = (I[97, :] * 0.3 + I[98, :] * 0.7 + I[100, :] * 0.7 + I[101, :] * 0.3) * 0.5
  32. I[219, :] = (I[217, :] * 0.3 + I[218, :] * 0.7 + I[220, :] * 0.7 + I[221, :] * 0.3) * 0.5
  33. 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
  34. 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
  35. # 对数操作
  36. I_nolog = I.copy()
  37. I0 = np.max(np.max(I))
  38. I = np.log(I / I0)
  39. I0 = np.max(np.max(I_noflat))
  40. I_noflat = np.log(I_noflat / I0)
  41. I0 = np.max(np.max(I_nobad))
  42. I_nobad = np.log(I_nobad / I0)
  43. # 反投影重建
  44. data = iradon(I, theta=list(range(360)), filter_name='ramp', interpolation='linear')
  45. noflat = iradon(I_noflat, theta=list(range(360)), filter_name='ramp', interpolation='linear')
  46. nobad = iradon(I_nobad, theta=list(range(360)), filter_name='ramp', interpolation='linear')
  47. nolog = iradon(I_nolog, theta=list(range(360)), filter_name='ramp', interpolation='linear')
  48. #显示结果
  49. plt.figure(2)
  50. plt.imshow(data, cmap='gray')
  51. plt.title('重建结果',fontproperties=myfont)
  52. plt.figure(3)
  53. plt.imshow(noflat, cmap='gray')
  54. plt.title('无平场校正',fontproperties=myfont)
  55. plt.figure(4)
  56. plt.imshow(nobad, cmap='gray')
  57. plt.title('无坏点补偿',fontproperties=myfont)
  58. plt.figure(5)
  59. plt.imshow(nolog, cmap='gray')
  60. plt.title('无对数操作',fontproperties=myfont)
  61. plt.show()

声明:本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:【wpsshop博客】
推荐阅读
相关标签
  

闽ICP备14008679号