CT骨料图像分割这件事，我从一开始的手动圈选到后来用Python脚本批量搞定，中间踩了不少坑。做混凝土细观模拟的人都懂，骨料的真实形状对力学性能影响很大，而CT扫描刚好能无损拿到颗粒的三维形态。问题在于，一坨试样扫下来少说几百张切片，单靠人工在二维图上一个一个描边，既不现实也没法保证再现性。后来我把整套流程搬到了Python上，用脚本串起“读图、降噪、分割、重建、导出STL”这条链路，千把个颗粒几分钟就能处理完，而且因为整个流程在本机跑，扫描数据不出实验室，稳定性和安全性都更有底。
写这篇东西，是想把这套方案完整讲清楚：为什么选Python这套技术栈、每一步的核心算法怎么选、批量导出时怎么保证效率和正确性。不管你是做混凝土、岩土还是其他多孔材料的CT图像分析，只要能拿到断层切片序列，这篇文章里的思路和代码基本都能直接拿过去改。
1. 处理思路与方案选型
1.1 为什么是Python + 本地批处理
早期我看到很多人处理CT骨料用的是商业软件，功能确实全，但在批量处理、结果可控和算法可定制这三个方面天然吃亏。商业软件通常是图形界面操作，对几十张切片还好，一旦换成几百个切片、上千个颗粒，手动点来点去容易出错，而且同样的数据过两周再做一次，参数没记录的话结果可能完全不一样，复现性是大问题。
Python这条路线的好处是流程可编程、参数可固化。numpy负责多维数组运算，scikit-image提供全套图像分割算法，trimesh处理三维网格和STL导出，这些库组合起来刚好覆盖从二维切片到三维模型的全部需求。更重要的是，整个脚本本地执行，不依赖云端服务，扫描数据不需要上传到任何平台，对于实验室数据来说，数据私密性这一块就踏实了。
我实际工作里最常见的诉求是：一批试样扫描完，想要逐颗粒的STL模型，给后续有限元网格划分或者3D打印用。这种批量提取需求，恰恰是脚本化处理的舒适区。写一次处理函数，循环跑一百个试样，中间不用人工干预，输出文件名和目录结构整齐划一，后面怎么查都方便。
1.2 CT切片到STL的完整数据链路
先梳理一下数据流。CT断层扫描得到的是一组二维灰度切片，沿着扫描轴向堆叠起来，就是三维体数据。这个体数据里，骨料因为密度高、CT值大，在灰度图上偏亮，水泥砂浆基体偏暗。我们要做的事，就是在三维体数据里把属于骨料的体素标出来，生成一个同样的三维掩膜（每个体素标记为0或1），然后从这块掩膜上提取表面网格，写成STL文件。
整个过程可以拆成五步，后面每一节都会展开讲：
- 读取切片序列，组装体数据；
- 体素尺寸标定，把像素坐标换成物理毫米坐标；
- 去噪与对比度增强，改善分割输入质量；
- 分割算法提取骨料掩膜，并分离粘连颗粒；
- 三维重建生成网格，批量导出STL。
很多人一上来就想用深度学习做分割，但实际应用中，规则算法往往更快更稳。骨料颗粒在CT图里和基体灰度差异明显，阈值加分水岭的组合拳，已经能解决绝大部分情况。深度模型参数多、训练成本高，对硬件有要求，对于这种形态相对规整、对比度足够的目标，我建议先用传统图像处理方法跑通，遇到难分的情况再考虑加模型不迟。
1.3 整体流程的五段式分解
回到具体的技术选型。在方案设计时我给自己定了几个原则：每段流程只做一件事，前一步的输出是后一步的输入；所有中间结果可以落盘，方便检查或断点续跑；算法参数放到配置区，不硬编码在代码里。
第一段是“数据组装”，把文件夹里几百张PNG或DICOM切片读成numpy三维数组。这一步别图省事用循环list拼接，内存峰值很高，最好先统计切片数量、宽高，预先分配numpy数组再填入数据。
第二段是“预处理”，中值滤波去噪、形态学操作清理孤立噪点。这里有个容易忽略的点：能不能在这一步就把造成分割困难的灰度不均匀问题解决掉。有些试样因为CT硬化伪影，中心区域灰度偏暗，边缘偏亮，直接全局阈值分割会整片误判，这时候需要先做偏置场校正或者局部归一化。
第三段是“分割”，先用Otsu自动求阈值得到料体和基体的二值图，再对二值图做连通域分析和距离变换，结合分水岭算法把相互接触的颗粒切开。
第四段是“后处理”，包括小体积噪点过滤、颗粒内部空腔填充、边缘毛刺清理。
第五段才是“重建与导出”，对每个颗粒掩膜单独跑marching cubes生成网格，再用网格简化和平滑控制面数，最后统一导出STL。导出的时候同时生成一个CSV清单，记录每个颗粒的编号、体积、表面积和等效球径，这个清单做统计分析特别好用。
2. 预处理：CT序列读取与像素尺度标定
2.1 从切片文件夹到numpy体数据
先解决把一堆二维切片变成三维数组的问题。我常用的方法是用imageio或cv2读取普通PNG序列，如果手头是DICOM格式，就搭配pydicom或者SimpleITK处理。这里给一段基础代码，演示从PNG序列组装体数据：
import numpy as np import imageio.v2 as iio import glob def load_slice_sequence(folder, sort_key='slice_'): files = sorted(glob.glob(folder + '/*.png'), key=lambda x: x.split(sort_key)[1]) # 先读一张确定形状 first = iio.imread(files[0]) h, w = first.shape[:2] volume = np.zeros((len(files), h, w), dtype=np.float32) for idx, f in enumerate(files): img = iio.imread(f) volume[idx] = img return volume注意这里有一个关键点：切片的排序不能只按文件名，得看扫描顺序。如果文件名里带编号但编号不是纯数字，字符串排序会出现“slice_10”排在“slice_9”前面的情况，所以排序逻辑上要提取出数字来比较。更稳妥的做法是用正则表达式从文件名里抽出序号再排序。
还有方向问题。CT切片从哪个方向开始扫描，体数据的z轴方向就对应哪个方向。如果后面要和有限元软件或者其他三维数据对齐，得搞清楚DICOM头里的方向余弦，否则模型在三维空间里可能是镜像或者翻转的。我一般的做法是先导出一个体数据和一个STL，在专业查看器里转一圈确认方向没问题，再批量跑后面的流程。
2.2 体素尺寸与单位标定
这一步是很多人容易忽略的，但它直接决定导出的STL模型物理尺寸对不对。CT切片里每个像素的实际大小由扫描分辨率决定，比如像素间距0.0732mm，切片间距0.5mm。如果你不把这些间距考虑进去，直接跑marching cubes，得到的模型在x和y方向是对的，但z方向会被拉长好几倍，整个颗粒形状完全变形。
切片间距（Slice Thickness）通常是DICOM头里的信息，PNG序列往往没有。这时候有两个选择：一是问扫描工程师要扫描参数，二是用已知尺寸的标定球做校准。我建议优先选用前者，省时省力。
组装体数据时，把这三个方向的间距记录成元组：
voxel_spacing_mm = (0.5, 0.0732, 0.0732) # (dz, dy, dx)这个元组要和numpy数组的shape顺序对应。volume数组的shape是(z, y, x)，所以spacing的顺序也必须是(z, y, x)。这个对应关系错一点，后面导出模型就会整体错位。
在重建阶段，把体素坐标缩放成物理坐标后，STL模型默认就是毫米单位。这部分代码在第四节里统一给出。
2.3 去噪与对比度归一化实操
CT图像里最常见的噪声是椒盐噪声和泊松噪声，还有因为射线硬化产生的条纹伪影。骨料分割对边缘的准确性要求高，所以滤波器的选择要同时照顾“去噪”和“保边”两个目标。中值滤波在去除孤立噪点的同时能较好地保留边缘，是性价比最高的一招，适合作为首选。
对三维体数据做中值滤波，直接拿scipy的medfilt写就行，但要注意窗口别太大。窗口太大会把细小的骨料颗粒整个磨掉，一般3x3x3或者5x5x5就够用了。切片之间如果有轻微的灰度不均匀，我还会在z方向单独做一层直方图匹配，把每一片的灰度分布拉齐，这会大大改善后续阈值分割的稳定性。
对比度归一化方面，如果整幅图的灰度范围特别宽，可以先做一个截断，比如把1%和99.9%的分位数之外的灰度裁掉，再做线性拉伸，让骨料和基体的对比更明显。这一步对阈值分割的影响很大，能有效减少阈值偏离带来的误差。
3. 分割算法：从灰度到颗粒掩膜
3.1 阈值分割的边界确定
骨料和砂浆基体的灰度差异是分割的基础。如果直方图呈现明显的双峰，用Otsu算法自动算出来的阈值基本可靠。Otsu原理不复杂——在所有灰度级别里找一个阈值，使得按它切分后两类的类内方差最小、类间方差最大。scikit-image里一行就能调：
from skimage.filters import threshold_otsu thresh = threshold_otsu(volume) # 全图统一阈值 mask = volume > thresh我个人经验是，先打印一下直方图，确认有没有明显双峰再跑Otsu。如果数据呈现单峰或者多峰交叠严重，Otsu会给出一个中间值，分割结果往往不太理想。这时候可以考虑多阈值方法，或者干脆手动指定一个范围，把灰度值介于[min_v, max_v]的体素判为骨料。
还有一点值得说：扫描时的电压电流参数会影响灰度绝对数值，不同批次的试样可能阈值不一样。所以脚本里阈值最好做成参数，跑之前先可视化几张检查一下，再批量执行。不要一套阈值跑到底，非常容易翻车。
3.2 粘连颗粒分离与分水岭算法
如果骨料颗粒相互独立，阈值分割后直接找连通域就够了。但混凝土试样里大量的骨料是互相挨着的，CT图上看过去一片亮区连着好几颗石子，直接连连通域会把它们当成一个大团块，导出STL之后形状和实际完全不符。
分离粘连颗粒的标准做法是距离变换加分水岭。先说原理：距离变换把二值图里的每个前景体素到最近背景的距离算出来，这样原本粘连的两个圆形颗粒，在距离图上会形成两个中心高、交界处低的“山包”。分水岭算法从这些山包的局部最高点开始“注水”，水漫到两个山包交界的低洼处时会形成分界线，也就是颗粒的分界面。
代码也不复杂：
from scipy import ndimage from skimage.feature import peak_local_max from skimage.segmentation import watershed distance = ndimage.distance_transform_edt(mask) coords = peak_local_max(distance, min_distance=5, exclude_border=2) markers = np.zeros_like(distance, dtype=int) markers[tuple(coords.T)] = np.arange(1, len(coords) + 1) labels = watershed(-distance, markers, mask=mask)这里面有两个参数要仔细调。min_distance控制的是两个种子点之间的最小距离，设置太小会把一个颗粒上的多个凸起误判成多个颗粒，设置太大又会把真正粘连的颗粒识别不了，合并成一个大块。exclude_border用于去掉边界上的假峰，避免试样边缘区域的噪点干扰。
watershed之前我把距离图取负，是因为分水岭默认找的是“山谷”而不是“山峰”，取负之后原来的山峰就变成了山谷，这样种子点就成了注水起点。
如果发现过分割（一个颗粒被切成好几块），可以用h-minima变换先对距离图做压制，过滤掉那些深度不够的局部极小值，再跑watershed。h值越大，过分割越少，但太大的话会把真实颗粒之间的浅凹槽也抹掉，颗粒又黏在一起。这个参数属于试出来最靠谱，我给的建议是先用小样本调，确认肉眼看起来合理再批量跑。
3.3 连通域分析与尺寸过滤
分离完成之后，labels里的每一个编号就是一个独立的颗粒。但这里面会混入一些噪点：要么是图像噪声产生的极小假目标，要么是试样切割边界上被截断的颗粒残片。后处理阶段可以通过体积过滤把它们清掉。
props = ndimage.measurements.label_stats(labels) for label_id in np.unique(labels): if label_id == 0: continue vol = int(np.sum(labels == label_id)) if vol * voxel_volume < min_vol_mm3: labels[labels == label_id] = 0其中voxel_volume是单个体素代表的体积，即z_spacing * x_spacing * y_spacing。min_vol_mm3按实际需求定，比如只关注粒径大于4.75mm的粗骨料，就可以按等效球径反算最小体积：4/3 * pi * (4.75/2)^3。颗粒粒径分布分析时，经常还要加一个最大长轴的过滤条件，这个用ndimage.measurements.center_of_mass和最小外接盒可以算，不过日常场景下体积过滤已经够用。
做完这些，掩膜质量基本就OK了。检查时可以在几个方向做最大密度投影，或者生成一个半透明三维渲染看一眼，确认没有大问题再进重建阶段。
4. 三维重建与STL批量导出
4.1 marching cubes提取表面网格
当我把掩膜数据弄好之后，重建这一步其实是最机械的。scikit-image的marching_cubes算法，可以把三维二值掩膜转换成三角面片网格。这个算法做的事情，通俗讲就是在体数据里找灰度值跨越等值面的位置，在这些位置插值生成顶点和三角面，把隐藏的表面“捞”出来。
关键参数是level。对于纯二值掩膜，取0.5通常能得到体积保持最好的结果。如果掩膜里已经是0和1，level=0.5意味着表面落在体素中心附近，不会过度膨胀也不会过度收缩。
from skimage.measure import marching_cubes verts, faces, normals, values = marching_cubes(mask.astype(np.uint8), level=0.5)拿到verts之后，现在要干的第二件事是把体素坐标转成物理毫米坐标。这一步不把第三步的spacing放进去的话，模型会变形得很厉害：
verts_mm = verts.copy() verts_mm[:, 0] *= voxel_spacing_mm[0] # z方向 verts_mm[:, 1] *= voxel_spacing_mm[1] # y方向 verts_mm[:, 2] *= voxel_spacing_mm[2] # x方向这里的顺序千万别搞混。marching_cubes返回的顶点坐标对应数组的(z, y, x)下标，所以和体数据的shape是一致的，必须和spacing元组一一对应。给每位颗粒生成网格时循环这个流程即可，每次处理一个颗粒的掩膜。
4.2 网格平滑与简化
直接生成的网格往往面片数量巨大，一个中等颗粒可能就有几十万个三角面。这样的STL文件动辄几十MB，对后续仿真软件读取和网格划分都不友好。两个后处理操作必不可少：平滑和简化。
平滑我用trimesh里的Laplacian平滑，迭代次数控制在10到20次。Laplacian平滑的原理是把每个顶点往邻居顶点的平均位置拉，多轮之后表面会变得细腻，但迭代太多会导致体积收缩，颗粒整体变小，细节被磨掉。我一般做完平滑前后各算一次mesh.volume，如果体积偏差超过2%，就说明迭代过量了，需要回调。
简化用的是二次误差测度（Quadric Error Metrics），trimesh里也有封装，通过mesh.simplify_quadric_decimation(target_faces)实现。这个算法的核心思想是每次合并边时，尽量让被移除顶点带来的几何误差最小，实际上就是反复合并一条边的两个端点，用二次型近似描述顶点附近曲面的变化。
import trimesh mesh = trimesh.Trimesh(vertices=verts_mm, faces=faces) mesh = mesh.simplify_quadric_decimation(target_faces=20000) mesh = mesh.laplacian_smooth(iterations=10, lamb=0.5)target_faces设多少，取决于用途。只想看看形状和统计体积，5000到10000面足够；要用来做有限元网格划分或者3D打印，建议保留20000到50000面，既保细节又控制文件体积。
平滑和简化的顺序也有讲究。先简化后平滑，可以在保证拓扑不崩的前提下把面数降下来，再统一磨掉棱角。如果先平滑后简化，会发现简化后表面又开始有微小褶皱，效果不如前者干净。
4.3 批量处理管线与文件命名规范
批量导出的核心是把上面的流程包装成一个函数，对每个颗粒掩膜调用一次，然后统一写文件。我这里给一段完整示例，展示颗粒级STL导出和统计清单生成：
import csv, os def export_particles(labels, spacing, out_dir, min_vol=10.0): os.makedirs(out_dir, exist_ok=True) rows = [] ids = np.unique(labels) ids = ids[ids != 0] voxel_vol = spacing[0] * spacing[1] * spacing[2] for i, pid in enumerate(ids): mask_i = (labels == pid) vol = mask_i.sum() * voxel_vol if vol < min_vol: continue verts, faces, _, _ = marching_cubes(mask_i.astype(np.uint8), level=0.5) verts[:, 0] *= spacing[0] verts[:, 1] *= spacing[1] verts[:, 2] *= spacing[2] mesh = trimesh.Trimesh(vertices=verts, faces=faces) mesh.fix_normals() # 可选：简化与平滑 if len(mesh.faces) > 50000: mesh = mesh.simplify_quadric_decimation(target_faces=20000) mesh = mesh.laplacian_smooth(iterations=10) mesh.export(os.path.join(out_dir, f"particle_{pid:04d}.stl")) rows.append([pid, vol, mesh.area, (6*vol/3.1415926)**(1/3)]) with open(os.path.join(out_dir, "summary.csv"), "w", newline="") as f: writer = csv.writer(f) writer.writerow(["id", "volume_mm3", "surface_mm2", "eq_diameter_mm"]) writer.writerows(rows)文件名我建议带上试样编号，比如sample01_agg_0001.stl，后面要追溯到哪个试样哪个颗粒，一眼就能看出来。summary.csv里记录体积、表面积和等效球径，做粒形统计分析时直接读这个文件，不用重新跑分割。
批量跑多个试样时，外层再加一个循环，遍历试样目录，调用load_slice_sequence和process_sample两个函数。这里有几个工程化要点：
- 每个试样处理之前先检查输出目录是否存在，已经处理过的试样就跳过，这样即使中途断了也不用从头再来；
- 所有中间结果（滤波后的体数据、分割掩膜）分别存成.npz或者.npy，方便随时回到任意一步调试；
- 日志记录每次跑批的参数和文件来源，结果出了问题还能追溯。
这套管线我实测跑过几百个切片的数据，单试样从读图到导出，一分钟上下能出几十个颗粒的STL，全程CPU占用稳定，内存峰值可控，属于真正能落地到日常工作流的方案。
5. 常见问题与排查技巧实录
5.1 颗粒粘连切不开
现象是导出的STL里出现“哑铃状”或者“葫芦状”的大块头，明明看着像两个颗粒却合并成了一个。我先检查分水岭参数，大概率是min_distance设太大，种子点数量不足，粘连区没有形成分界线。解决步骤是：先跑一次peak_local_max，把找到的种子点数打印出来，和目视判断的颗粒数量对比，如果明显偏少就把min_distance调小2到3个像素，重新跑。
还有一种情况是阈值分割阶段就把两个颗粒连成一片了，中间只有很窄的一点亮的连接。这时候光调分水岭效果有限，建议在阈值分割后先对mask做一次开运算（先腐蚀后膨胀），把细小的连接通道断开，再跑距离变换和分水岭。
另外，颗粒表面粗糙、有棱角时，距离变换在颗粒内部可能出现很多局部高点，导致一个颗粒被切碎。这种情况用h-minima压制或提高min_distance都能改善，我一般两个参数一起调。
5.2 STL表面孔洞与碎片
导出网格后拿查看器打开，发现表面有孔洞、或者有悬浮的小三角形碎片飞在外面。孔洞通常是掩膜表面本身有薄弱区域，或者marching_cubes在低对比度区域插值不完整。处理的方式是用trimesh的fill_holes方法，把网格的边界环补成三角面。
碎片问题更多是掩膜后处理没做干净。我插过一两次手，总结出来的办法是：在连通域分析时，把每个候选颗粒的物理体积算出来，低于最小体积阈值的直接剔除，不给它生成网格的机会。还有一个隐蔽问题：如果掩膜边界贴在图像边缘，切出来的颗粒在边界外是敞口的，marching_cubes生成的网格不是闭合的，导出后一看是个“碗”。处理方法是给体数据四周加一圈背景padding，再做连通域和重建，最后再把偏移坐标减去。
5.3 模型尺寸与物理尺寸对不上
这个问题出现的频率比想象中高。要么是切片间距填错，要么是单位约定不一致，要么是方向搞反了。我的排查顺序是：先在体数据里量一个已知颗粒的体素数，乘上对应的spacing，和实际测量值对比，确认单一方向对不对；再导出一个简单形状（最好是球体或者圆柱体）的STL，在查看器里测尺寸，确认整体缩放。
第二个容易踩的坑是STL只有顶点坐标、没有单位信息，不同软件默认单位可能不一样。常见约定是STL文件用毫米单位，如果后续要导入到以米为单位的仿真软件，需要换算，这个换算在导出时做掉，不要等下游软件再去改。
5.4 内存占用过高与处理速度慢
体积数据动辄几百MB，如果全程用float64存，内存很容易爆。优化手段有三条：一是读图时就用float32甚至uint16存储，二是在每个阶段处理完后释放中间变量，三是处理时可以分块。CT骨料分割的主要计算瓶颈在距离变换和marching cubes，这两个操作对内存都是“全宽度”的，最好确保机器至少有16GB内存，大试样能到32GB更好。
如果数据实在太大，还有一个思路是把三维分割拆成“粗分割+局部精化”。先用降采样的体数据快速找到每个颗粒的大致区域，再在每个颗粒的局部邻域里跑高精度分割和重建。这样处理时间能大幅下降，内存占用也从“整个体数据”降到“单颗粒局部邻域”。
5.5 骨料与基体灰度接近时的分割难题
最后补充一个偏研究向的问题。有些骨料（尤其是石灰岩类）和水泥基体在CT图里的灰度差异很小，全局阈值分割经常会丢失部分颗粒或者把基体误判成骨料。常规手段是一遍一遍调阈值，但效果有限。
我做过一次之后的心得是：遇到这种情况，优先考虑两步分割。第一步先用较宽松的阈值把“确定是骨料”的核心区提出来，第二步在剩余区域用边缘检测或者局部自适应阈值补齐边界。形态学的闭运算也能帮忙，把颗粒内部的孔隙填上，让掩膜更完整。实在还不行，再考虑训练一个轻量的分类模型，输入换成多尺度灰度特征，输出逐体素的骨料概率。这一步能用传统方法解决就不上模型，逻辑更简单，结果也更容易解释。
最后说点实际操作体会
整套流程写下来，看起来都是现成算法，但真正跑通一个可靠的生产管线，要在很多细节上花时间。我最早自己实现的时候，直接在完整体数据上跑marching cubes，结果一个颗粒没分出来，找半天才发现是阈值分割阶段就把整片基体误判成了骨料，滤波参数还没调对。后来学乖了，每一步都保存中间结果，跑几步停下来用专业查看器检查，确认没问题再继续。
调参上我有一个习惯，先用二三十张切片组成的“最小可验证子集”跑通全流程，这个子集不用完整覆盖试样，只要包含典型情况就行。全流程跑通、参数确定之后，再一次性放全量数据下去跑。这样排错成本低很多，也避免一遍遍重复跑全量遮挡问题。
后面如果要扩展，这个流程还能加很多功能：按颗粒粒径分布自动分类统计、富集区域筛选、多阈值联合分割、甚至结合传统分割和简单模型提升难例效果。至少目前这套Python方案，已经能稳定地把批量CT数据的骨料分割和STL导出做成一个自动化、可复现的流程。希望你看了之后跑自己的数据能少踩几个坑，拿着代码跑通一遍，思路就全通了。
