## 1. 从零开始：认识ASTRA Toolbox与Python断层扫描 如果你对医学影像、工业CT或者材料科学里的三维成像感兴趣，那你很可能听说过“断层扫描重建”这个词。简单来说，这就像给一个物体拍很多张不同角度的“X光片”，然后通过计算机算法，把这些二维的“影子”拼凑还原成物体内部的三维结构。这个过程听起来很酷，但实际操作起来，尤其是自己动手写代码实现，往往会遇到一堆让人头疼的数学和编程问题。 几年前，我在处理一批生物样本的微CT数据时就深有体会。当时试了好几个开源工具，要么速度慢得让人抓狂，要么配置复杂得像在解谜。直到我遇到了 **ASTRA Toolbox**，才感觉终于找到了趁手的“兵器”。ASTRA Toolbox是一个专门为X射线断层扫描（包括平行束、锥束等）设计的开源算法工具箱，它的核心优势在于**高效**和**灵活**。高效，是因为它底层大量使用了GPU（CUDA）进行加速，处理速度比纯CPU实现快几十倍甚至上百倍，这对于动辄几百上千张投影图像的重建来说，简直是救命稻草。灵活，是因为它提供了从底层几何定义到高层算法调用的完整Python接口，让你不仅能“用”，还能“改”和“创”，真正理解重建过程的每一个细节。 那么，谁适合学习使用它呢？我觉得主要有三类朋友：一是**科研人员**，比如做生物医学成像、材料无损检测的，需要快速验证算法或处理实验数据；二是**工程师**，从事工业CT系统开发或图像处理算法优化；三是**有一定Python基础的学生或爱好者**，想深入理解计算机断层扫描的原理并动手实践。如果你属于其中任何一类，那么这篇实战指南就是为你准备的。我会把我踩过的坑、总结的经验，用最直白的方式分享出来，让你能快速上手，用Python和ASTRA Toolbox玩转X射线重建。 ## 2. 环境搭建与核心概念初探 工欲善其事，必先利其器。在写第一行重建代码之前，我们需要先把环境准备好。ASTRA Toolbox的安装不算特别复杂，但有几个关键点需要注意，否则很容易卡在第一步。 ### 2.1 安装与配置：避开我踩过的那些坑 ASTRA Toolbox的安装主要依赖两个部分：一是工具箱本身的Python包，二是用于GPU加速的CUDA环境。官方推荐使用`conda`来管理环境，这确实能省去很多麻烦。 首先，创建一个新的conda环境是个好习惯，能避免包版本冲突。打开你的终端（或Anaconda Prompt），执行： ```bash conda create -n astra_env python=3.8 conda activate astra_env ``` 这里我选择Python 3.8，是因为它在与一些科学计算库的兼容性上比较稳定。当然，3.9或3.10通常也可以。 接下来是安装ASTRA Toolbox。最直接的方式是通过conda-forge频道安装： ```bash conda install -c astra-toolbox/label/dev -c conda-forge astra-toolbox ``` 注意，我们这里指定了`label/dev`频道，这能确保安装到包含最新功能和修复的版本。安装过程会自动处理许多依赖，比如NumPy、SciPy等。 **但是，最重要的部分来了：GPU支持。** ASTRA的强大性能很大程度上依赖于CUDA。如果你有一张NVIDIA显卡并希望使用GPU加速（强烈建议），你必须确保系统上安装了正确版本的CUDA工具包。ASTRA Toolbox通常与特定版本的CUDA绑定。例如，当前版本可能兼容CUDA 11.x。你可以在安装时指定CUDA版本： ```bash conda install -c astra-toolbox/label/dev -c conda-forge astra-toolbox cudatoolkit=11.2 ``` 安装完成后，强烈建议运行一个简单的测试来验证安装是否成功，以及GPU是否可用： ```python import astra print(astra.astra.use_cuda()) # 如果输出True，恭喜，GPU可用！ ``` 如果这里输出`False`，别慌。可能是CUDA路径问题或者驱动不匹配。你可以尝试在Python中检查`astra.astra.version()`，或者去ASTRA的GitHub仓库的issue区看看，很多常见问题都有解决方案。我当初就花了半天时间折腾CUDA和驱动版本，所以这一步耐心点很值得。 ### 2.2 理解两个核心：几何与数据对象 环境搞定后，我们先别急着写重建代码。要想用好ASTRA，必须理解它的两个核心概念：**几何（Geometry）**和**数据对象（Data Object）**。你可以把整个重建过程想象成拍电影：几何定义了“摄影棚”的布局（相机在哪，怎么移动），数据对象则是“胶片”和“最终成片”。 **几何（Geometry）** 分为两大类： 1. **体积几何（Volume Geometry）**：它定义了你要重建的那个“物体”所在的空间盒子。这个盒子是三维的（对于3D重建），你只需要告诉ASTRA这个盒子在Y, X, Z三个方向上有多少个“小格子”（体素，voxel）。例如，`vol_geom = astra.create_vol_geom(256, 256, 256)`就定义了一个边长为256体素的正方体空间，你的目标就是求出这个正方体内每一个小格子的密度值。 2. **投影几何（Projection Geometry）**：它定义了X射线源和探测器是如何围绕物体运动的。这是重建中最容易出错的部分！ASTRA支持多种几何，比如最常用的锥束（Cone beam）。创建一个锥束几何需要一堆参数： ```python proj_geom = astra.create_proj_geom( 'cone', # 几何类型：锥束 1.0, 1.0, # 探测器单个像素的宽和高（单位自定，但需一致） 300, 400, # 探测器面板的行数和列数（即探测器的分辨率） angles, # 一个数组，定义了每个投影角度（弧度制） 500, # 射线源到物体旋转中心的距离 300 # 物体旋转中心到探测器平面的距离 ) ``` 这里的`angles`通常用`np.linspace(0, 2*np.pi, 360, endpoint=False)`生成，表示在0到360度范围内均匀采集360个角度的投影。`source_origin`和`origin_det`这两个距离参数至关重要，它们决定了投影的放大倍数和几何畸变，一定要根据你的实际实验设备或仿真设定来设置。 **数据对象（Data Object）** 是ASTRA内部管理数据的方式。我们不会直接操作NumPy数组，而是先创建数据对象，把数组“喂”给它，然后通过一个ID来引用它。这听起来有点绕，但这样做是为了高效地在CPU和GPU内存之间传输数据。主要用到的函数是`astra.data3d.create`： - `sino_id = astra.data3d.create('-proj3d', proj_geom, sinogram_data)`：创建一个投影数据对象。`sinogram_data`是一个三维NumPy数组，形状通常是`(探测器行数, 投影角度数, 探测器列数)`。 - `vol_id = astra.data3d.create('-vol', vol_geom, volume_data)`：创建一个体积数据对象。`volume_data`是一个三维NumPy数组，形状与`vol_geom`定义的尺寸一致。 创建对象后，你可以用`astra.data3d.get(object_id)`把数据取回成NumPy数组查看。**切记**，所有创建的对象在用完后，要用`astra.data3d.delete(object_id)`手动删除，尤其是在循环或大量数据处理中，否则会导致内存泄漏。 ## 3. 实战第一步：从模拟物体生成投影数据 在真正处理实验数据之前，用一个已知的模拟物体（比如一个立方体、一个球）来生成投影数据，是验证整个流程是否正确的黄金标准。这就像在调试电路时，先输入一个标准信号看看输出对不对。 ### 3.1 构建一个“数字幻影” 我们先来造一个简单的三维模型——一个悬浮在空中的立方体。在断层扫描领域，这被称为“数字幻影”（Digital Phantom）。 ```python import numpy as np import astra # 1. 定义体积几何：一个128x128x128的空间 vol_geom = astra.create_vol_geom(128, 128, 128) # 2. 创建幻影数据：初始化为全0（代表空气或背景） phantom = np.zeros((128, 128, 128), dtype=np.float32) # 3. 在中心区域画一个立方体：将30:100范围内的体素值设为1（代表高密度材料） phantom[30:100, 30:100, 30:100] = 1.0 # 可视化一下中间切片看看 import matplotlib.pyplot as plt plt.figure(figsize=(6,6)) plt.imshow(phantom[64], cmap='gray') # 取Z轴中间的切片 plt.title('Digital Phantom (Central Slice)') plt.colorbar() plt.show() ``` 运行这段代码，你会看到一个正方形的白色方块出现在灰色背景中央。这个“1.0”的值代表的是相对衰减系数，在实际物理中对应某种材料的X射线吸收能力。 ### 3.2 配置几何并执行正投影 有了“物体”，现在我们来布置“拍摄现场”。我们模拟一个锥束CT系统，射线源和探测器绕物体旋转一周。 ```python # 1. 定义投影角度：360度内均匀取180个角度（通常180度理论上就够了，但这里为了演示用360度） angles = np.linspace(0, 2 * np.pi, 180, endpoint=False) # 2. 定义探测器参数 det_row_count = 256 # 探测器行数（纵向像素） det_col_count = 256 # 探测器列数（横向像素） det_spacing_x = 1.0 # 探测器像素宽度 det_spacing_y = 1.0 # 探测器像素高度 # 3. 定义系统几何参数（这些值需要根据实际或仿真设定调整） source_origin = 300 # 射线源到旋转中心的距离 origin_det = 300 # 旋转中心到探测器的距离 # 4. 创建锥束投影几何 proj_geom = astra.create_proj_geom('cone', det_spacing_x, det_spacing_y, det_row_count, det_col_count, angles, source_origin, origin_det) ``` 接下来，就是激动人心的“拍摄”环节——正投影（Forward Projection）。这个过程在ASTRA里可以通过创建一个“正投影算法”对象来完成。 ```python # 1. 创建投影数据对象（此时内容为空） proj_id = astra.data3d.create('-proj3d', proj_geom) # 2. 创建体积数据对象，并将我们的幻影数据载入 vol_id = astra.data3d.create('-vol', vol_geom, data=phantom) # 3. 配置正投影算法（使用GPU加速的FP3D_CUDA） cfg = astra.astra_dict('FP3D_CUDA') cfg['VolumeDataId'] = vol_id # 指定输入物体 cfg['ProjectionDataId'] = proj_id # 指定输出投影数据存放处 proj_alg_id = astra.algorithm.create(cfg) # 4. 执行算法，生成投影数据！ print("开始正投影计算...") astra.algorithm.run(proj_alg_id) print("正投影计算完成！") # 5. 获取生成的投影数据（sinogram） sinogram = astra.data3d.get(proj_id) print(f"投影数据形状: {sinogram.shape}") # 应该是 (256, 180, 256) # 6. 可视化第90个角度的投影 plt.figure(figsize=(12,5)) plt.subplot(121) plt.imshow(phantom[64], cmap='gray') plt.title('Original Phantom (Slice)') plt.subplot(122) plt.imshow(sinogram[90, :, :], cmap='gray') # 查看第90个角度的投影图像 plt.title('Projection at Angle 90') plt.colorbar() plt.show() ``` 当你看到右边那张由许多平行线组成的投影图时，就说明正投影成功啦！这张图就是你的立方体在某个特定角度下的“X光片”。`sinogram`这个三维数组里，存储了物体在所有180个角度下的投影。 ### 3.3 数据保存与格式转换 生成的数据最好保存下来，方便后续重建或分析。我推荐使用`.npy`格式保存NumPy数组，它比图像序列更紧凑，且能保留完整的浮点数精度。 ```python # 保存投影数据 np.save('my_sinogram.npy', sinogram) # 同时保存几何参数（因为重建时需要完全相同的几何！） import pickle geometry_info = { 'proj_geom': proj_geom, 'vol_geom': vol_geom, 'angles': angles } with open('geometry_params.pkl', 'wb') as f: pickle.dump(geometry_info, f) ``` **最后，也是最重要的好习惯：清理现场。** ASTRA创建的所有对象ID都会占用内存（尤其是GPU内存），必须手动删除。 ```python astra.algorithm.delete(proj_alg_id) astra.data3d.delete(proj_id) astra.data3d.delete(vol_id) ``` ## 4. 核心环节：使用不同算法进行图像重建 投影数据到手后，就进入了最核心的图像重建环节。ASTRA提供了从解析法到迭代法的多种算法，各有优劣，适用于不同场景。 ### 4.1 快速解析法：FDK算法上手 对于锥束几何，最经典、最快的重建算法是FDK算法（Feldkamp-Davis-Kress）。它是一种滤波反投影算法，速度极快，通常作为第一选择。 ```python # 假设我们已经从文件加载了投影数据 sinogram 和几何参数 proj_geom, vol_geom # 1. 创建用于存放重建结果的数据对象 recon_id = astra.data3d.create('-vol', vol_geom) # 初始内容为0 # 2. 配置FDK算法（GPU版） cfg_fdk = astra.astra_dict('FDK_CUDA') # 注意算法名是'FDK_CUDA' cfg_fdk['ProjectionDataId'] = proj_id cfg_fdk['ReconstructionDataId'] = recon_id # FDK算法通常不需要太多额外参数，配置非常简单 # 3. 创建并运行算法 fdk_alg_id = astra.algorithm.create(cfg_fdk) astra.algorithm.run(fdk_alg_id) # FDK是解析法，一次执行即可完成重建 # 4. 获取重建结果 reconstruction_fdk = astra.data3d.get(recon_id) # 5. 可视化重建结果的一个切片 plt.figure(figsize=(8,4)) plt.subplot(121) plt.imshow(phantom[64], cmap='gray') plt.title('Original Phantom') plt.subplot(122) plt.imshow(reconstruction_fdk[64], cmap='gray') plt.title('FDK Reconstruction') plt.colorbar() plt.tight_layout() plt.show() ``` FDK重建速度非常快，对于我们的模拟数据，可能一瞬间就完成了。对比左右两图，你应该能看到重建出的立方体。但是，你可能会发现重建的图像边缘没有原图那么锐利，或者背景有少许噪声和伪影（比如条纹）。这是FDK算法的特点：**快，但在投影数据不完整、有噪声或者几何不理想时，重建质量会下降。** ### 4.2 高质量迭代法：SIRT算法深度解析 当数据质量不高（比如低剂量、少角度、有噪声）时，迭代重建算法往往能获得更好的结果。其中，SIRT（同步迭代重建技术）是一个在稳定性和效果上取得很好平衡的算法。 ```python # 1. 同样，创建重建数据对象 recon_id_sirt = astra.data3d.create('-vol', vol_geom) # 2. 配置SIRT算法（GPU版） cfg_sirt = astra.astra_dict('SIRT3D_CUDA') cfg_sirt['ProjectionDataId'] = proj_id cfg_sirt['ReconstructionDataId'] = recon_id_sirt # 可以设置迭代次数和约束条件 cfg_sirt['option'] = { 'MinConstraint': 0.0, # 体素最小值约束（物理上密度非负） 'MaxConstraint': 1.5 # 体素最大值约束（根据先验知识设定） } # 3. 创建算法对象 sirt_alg_id = astra.algorithm.create(cfg_sirt) # 4. 运行迭代。迭代次数是SIRT的关键参数！ num_iterations = 50 print(f"开始SIRT迭代重建，共{num_iterations}次迭代...") astra.algorithm.run(sirt_alg_id, iterations=num_iterations) print("SIRT重建完成！") # 5. 获取结果 reconstruction_sirt = astra.data3d.get(recon_id_sirt) ``` SIRT算法需要迭代多次，每次迭代都会根据当前重建图像与原始投影数据的差异进行更新。迭代次数太少，重建不充分；太多，则可能过拟合噪声且计算时间长。**如何选择迭代次数？** 一个实用的方法是观察重建误差随迭代次数的变化。我们可以稍微修改代码，在每次迭代后记录误差： ```python # 创建一个列表记录每次迭代后的数据误差（这里用投影数据误差的范数近似） errors = [] for i in range(1, 101): astra.algorithm.run(sirt_alg_id, 1) # 每次运行1次迭代 # 获取当前重建结果 recon_current = astra.data3d.get(recon_id_sirt) # (此处简化，实际计算误差需要重新正投影并与原始投影比较，略复杂) # 我们可以简单计算重建图像的变化率作为参考 if i > 1: diff = np.linalg.norm(recon_current - recon_previous) errors.append(diff) recon_previous = recon_current.copy() plt.plot(range(2, 101), errors) plt.xlabel('Iteration Number') plt.ylabel('Change per Iteration') plt.title('SIRT Convergence Behavior') plt.grid(True) plt.show() ``` 当曲线变得平缓，说明算法已经收敛，再增加迭代次数收益不大。对于我们的模拟数据，可能50-100次迭代就够了。对比FDK和SIRT的结果，你会发现SIRT重建的图像通常更平滑，噪声和伪影更少，但计算时间也长得多。 ### 4.3 算法对比与参数调优指南 为了更直观地对比，我们可以把FDK、SIRT（20次迭代）、SIRT（100次迭代）的结果放在一起看。 ```python fig, axes = plt.subplots(2, 2, figsize=(10, 10)) im0 = axes[0,0].imshow(phantom[64], cmap='gray') axes[0,0].set_title('Ground Truth') plt.colorbar(im0, ax=axes[0,0]) im1 = axes[0,1].imshow(reconstruction_fdk[64], cmap='gray') axes[0,1].set_title('FDK Reconstruction') plt.colorbar(im1, ax=axes[0,1]) im2 = axes[1,0].imshow(reconstruction_sirt_20[64], cmap='gray') # 假设有20次迭代的结果 axes[1,0].set_title('SIRT (20 iterations)') plt.colorbar(im2, ax=axes[1,0]) im3 = axes[1,1].imshow(reconstruction_sirt_100[64], cmap='gray') # 假设有100次迭代的结果 axes[1,1].set_title('SIRT (100 iterations)') plt.colorbar(im3, ax=axes[1,1]) plt.tight_layout() plt.show() ``` | 算法 | 优点 | 缺点 | 适用场景 | | :--- | :--- | :--- | :--- | | **FDK_CUDA** | **速度极快**，单次计算即可完成；实现简单，参数少。 | 对数据完备性要求高；噪声放大明显；锥束伪影可能较重。 | 数据质量好（高信噪比、全角度）、需要快速预览结果的场景。 | | **SIRT3D_CUDA** | **抗噪声能力强**；重建质量稳定，伪影少；可加入物理约束（如非负性）。 | **速度慢**，需要多次迭代；迭代次数需要调优；内存占用相对较高。 | 低剂量CT、稀疏角度CT、数据有噪声或缺失等不理想情况。 | | **CGLS3D_CUDA** | 收敛速度通常比SIRT快；对于某些问题更高效。 | 对条件数敏感，可能不稳定；参数（如预处理子）设置更复杂。 | 需要快速迭代收敛且问题条件较好的情况。 | | **FBP3D_CUDA** | 平行束几何下的标准解析法，理论完备。 | **仅适用于平行束几何**，锥束需用FDK。 | 传统的平行束CT扫描数据。 | **参数调优小贴士**： - **FDK**：关键参数是重建时的滤波器。在`cfg['option']`中可以设置`FilterType`（如‘Ram-Lak‘, ‘Shepp-Logan‘, ‘Cosine‘）和`FilterParameter`（截止频率）。`Ram-Lak`滤波器保留高频细节多，但噪声也大；`Shepp-Logan`或`Cosine`更平滑。 ```python cfg_fdk['option'] = {'FilterType': 'Shepp-Logan', 'FilterParameter': 0.8} ``` - **SIRT**：除了迭代次数，`MinConstraint`和`MaxConstraint`非常有用。如果你知道被测物体的密度范围（比如在0到某个值之间），设置约束可以显著改善重建质量，防止出现不合理的负值或过高的值。 - **通用技巧**：如果重建结果出现严重的环形伪影，可能是探测器像素增益不一致或几何标定不准。可以尝试在投影数据预处理阶段进行平场校正（Flat-field correction）。如果图像模糊，可以尝试在迭代算法中引入**总变分（TV）正则化**，ASTRA也提供了相关的算法，能更好地保留边缘同时抑制噪声。 ## 5. 处理真实数据与高级技巧 掌握了模拟数据的流程后，我们就可以挑战真实数据了。真实数据往往伴随着噪声、畸变和不完备性，这就需要一些额外的处理技巧。 ### 5.1 导入与预处理真实投影数据 真实数据通常来自CT扫描仪，保存为一系列TIFF或DICOM图像，或者一个专用的数据文件。我们需要将它们读入并整理成ASTRA需要的三维数组格式 `(探测器行数, 投影角度数, 探测器列数)`。 ```python import tifffile as tiff import os # 假设投影图像按顺序保存在‘projections/‘文件夹下，文件名为‘proj_001.tif‘, ‘proj_002.tif‘... data_dir = 'projections/' file_list = sorted([f for f in os.listdir(data_dir) if f.endswith('.tif')]) num_angles = len(file_list) # 读取第一张图像获取尺寸 first_img = tiff.imread(os.path.join(data_dir, file_list[0])) det_row_count, det_col_count = first_img.shape # 初始化数组 sinogram_real = np.zeros((det_row_count, num_angles, det_col_count), dtype=np.float32) # 循环读取所有图像 for i, filename in enumerate(file_list): img = tiff.imread(os.path.join(data_dir, filename)) sinogram_real[:, i, :] = img.astype(np.float32) print(f"真实投影数据形状: {sinogram_real.shape}") ``` 读入数据后，**预处理至关重要**。常见的步骤包括： 1. **对数变换**：CT探测器测量的是穿透后的光强I，而重建需要的是衰减系数 μ，其关系为 μ ∝ -log(I/I0)。其中I0是空白扫描（Flat-field）的光强。 ```python # 假设我们已经有了空白扫描图像 flat_field 和暗场图像 dark_field # 通常每个角度都有对应的空白和暗场，这里简化为全局校正 I0 = flat_field - dark_field I = sinogram_real - dark_field # 防止除零或log(0) I0[I0 <= 0] = 1e-6 I[I <= 0] = 1e-6 sinogram_corrected = -np.log(I / I0) ``` 2. **坏点校正**：探测器可能有死像素或响应异常的像素，可以用邻近像素插值替换。 3. **中心偏移校正**：旋转轴的中心必须精确对准，否则重建图像会模糊。ASTRA的投影几何中，默认旋转中心在体积几何的中心。如果实际数据有偏移，需要在创建`proj_geom`时，通过额外的参数（如`det_center_x`）进行设置，或者对投影数据进行平移。 ### 5.2 应对挑战：少角度与低剂量重建 真实实验中，为了减少扫描时间或辐射剂量，我们常常面临“少角度”或“低剂量”问题。投影数据不足，FDK这类解析法就会产生严重的条纹伪影。这时，迭代法结合**正则化**技术是更优的选择。 ASTRA内置了基于**总变分（TV）** 正则化的算法，如`SIRT-TV`。TV正则化倾向于寻找一个“分段常数”的解，即图像由大片均匀区域和清晰的边界构成，这非常符合很多实际物体的特性。 ```python # 配置SIRT-TV算法（可能需要从插件或特定版本获取，具体函数名请查文档） # 以下配置示意其思路 cfg_sirt_tv = astra.astra_dict('SIRT3D_CUDA_TV') # 注意算法名可能不同 cfg_sirt_tv['ProjectionDataId'] = proj_id cfg_sirt_tv['ReconstructionDataId'] = recon_id cfg_sirt_tv['option'] = { 'MinConstraint': 0.0, 'TVLambda': 0.01, # TV正则化的强度参数，需要调试 'TVIterations': 5 # 内部TV子问题的迭代次数 } ``` `TVLambda`参数控制平滑强度：值太大，图像会过于平滑，细节丢失；值太小，去噪和抑制伪影的效果不明显。这通常需要通过试错，或者基于一些图像质量指标（如信噪比、结构相似性）来调整。 ### 5.3 结果可视化与定量分析 重建出三维体数据后，我们不仅想看一个切片，还想从各个角度观察，甚至进行定量分析。 ```python # 1. 多切片查看 fig, axes = plt.subplots(2, 3, figsize=(15, 10)) slice_indices = [30, 64, 98] # 查看不同深度的切片 for i, idx in enumerate(slice_indices): axes[0, i].imshow(reconstruction_sirt[idx], cmap='gray') axes[0, i].set_title(f'Slice Z={idx}') axes[0, i].axis('off') axes[1, i].hist(reconstruction_sirt[idx].flatten(), bins=50) axes[1, i].set_title(f'Histogram at Z={idx}') axes[1, i].set_xlabel('Voxel Value') axes[1, i].set_ylabel('Count') plt.tight_layout() plt.show() # 2. 三维等值面渲染（使用matplotlib的3D功能或Mayavi、PyVista等库） # 这里使用matplotlib的简单3D可视化，适合小数据预览 from mpl_toolkits.mplot3d import Axes3D # 创建一个阈值，只显示密度大于0.5的区域 threshold = 0.5 z, y, x = np.where(reconstruction_sirt > threshold) fig = plt.figure(figsize=(10,8)) ax = fig.add_subplot(111, projection='3d') ax.scatter(x, y, z, c=reconstruction_sirt[z, y, x], cmap='hot', s=1, alpha=0.6) ax.set_xlabel('X') ax.set_ylabel('Y') ax.set_zlabel('Z') ax.set_title('3D Scatter Plot of High-Density Regions') plt.show() ``` 直方图可以帮助你分析重建值的分布，判断背景和材料的分离程度。三维散点图能给你一个物体空间结构的直观感受。对于更专业的三维渲染，我推荐使用**PyVista**或**Napari**库，它们能交互式地浏览体数据，功能强大得多。 处理真实数据时，我最大的体会是：**预处理和几何标定的重要性，不亚于重建算法本身。** 一个小的几何参数误差，就可能导致整个重建失败。因此，在扫描实物前，用已知的标准样品（如一个均匀小球）进行系统标定，是保证结果可靠的关键一步。ASTRA Toolbox给了我们强大的算法武器，但要打出精准的“十环”，还需要我们对整个物理成像过程有细致的理解和把控。