在实验室里处理粉末、颗粒类样品的人几乎都经历过同一件事坐在显微镜前一张图一张图地手动数颗粒、量尺寸。颗粒多了眼睛发花颗粒不规则了又不知道按哪个方向量更别说不同人数出来的结果能差出一截。后来我搭了一套基于 MATLAB 的不规则颗粒粒径、周长、面积测量及计数系统把“读图—分割—计数—量测—汇总”整条链路做成自动化流程一批样品几分钟就能出一份统计表。这篇就完整拆解这套系统的设计思路、核心算法、完整代码以及我这几年反复调参踩出来的坑。这套东西适合谁看主要是需要处理粉末冶金、矿物加工、陶瓷原料、制药颗粒、土壤或微塑料这类不规则颗粒图像的人。不管你是把颗粒拍成了显微照片还是用体式镜、扫描电镜出图只要图像里颗粒和背景能分得开这套处理框架基本都能用。不用懂太多图像处理理论跟着每个步骤把代码跑通就能得到一组可复现、可追溯的粒径统计结果。1. 先别急着写代码你要测的“粒径”到底是哪一个1.1 不规则颗粒为什么不能套用规则几何公式我见过不少人拿到颗粒图第一反应是“量直径”。等到颗粒形状不规则的时候这个问题就变得很尴尬一个长条形的颗粒对角线量是 80 个像素宽度量只有 30 个像素你说它的“直径”是多少规则物体能用一个半径或直径概括是因为它的形状信息可以由一个几何参数完全决定。但不规则颗粒的边界是任意曲线不存在这个性质。更关键的是图像上得到的只是颗粒的二维投影不是三维实体。也就是说你根本没用量到体积只能量到投影面积和投影周长。基于这两项再通过约定俗称的换算关系推算出各种“等效粒径”。所以在做系统之前先要承认一个事实图像测量给不出一个唯一的“真粒径”它给的是在特定定义下的等效值。这才是一切后续统计的地基。1.2 主流的几种等效粒径定义我把工程和科研里常见的粒径定义整理了一张表它们是同一颗颗粒在不同口径下的不同数字指标定义典型用途对颗粒摆放方向的敏感性等效圆直径ECD与颗粒投影面积相等的圆直径D√(4A/π)通用粒径统计、粒度分布不敏感费雷特径Feret沿某个方向作两条平行切线两线之间的最大/最小/平均间距筛分分析、取向分析敏感马丁径Martin沿固定方向将颗粒投影面积等分的弦长纤维、棒状颗粒分析敏感等效椭圆长短轴与颗粒二阶矩等效的椭圆的长轴、短轴形状与长径比评价不敏感圆形度C4πA/P²圆为 1越不规则数值越小颗粒棱角、破碎程度评价不敏感同一颗长度 80 像素、宽度 30 像素的长条形颗粒它的最大费雷特径可能接近 80 像素而等效圆直径只有大约 40 像素。这两个数字并没有谁对谁错只是回答的问题不一样。做筛分对比的人更关心费雷特径因为筛孔卡的是颗粒的宽度方向做比表面积或沉降分析的人更关心等效圆直径。1.3 这套系统最终采纳的指标口径我最终在系统里输出的核心指标是投影面积、等效圆直径、周长、长轴与短轴、圆形度另外附上最大费雷特径作为辅助参考。选择这个组合的理由很实际等效圆直径对颗粒的摆放方向不敏感同一颗颗粒在视野里转任意角度测出来的值都稳定适合做批间对比。周长本身单独用价值不大但它是计算圆形度所必需的。圆形度直接反映颗粒破碎、磨损的程度在选矿、研磨工艺里是特别关键的指标。长轴与短轴用于评估颗粒是否呈针状或片状这个信息在制药结晶、陶瓷粉体领域非常常用。需要额外说明的是这套系统测的不是严格意义上的三维粒度。假如你的样品是球形颗粒二维投影直径可以近似等于三维直径但如果是片状或针状颗粒投影测量会系统性地偏离真实三维尺寸这点需要在报告里明确标注否则下游工艺人员很容易误读数据。2. 图像预处理决定这套系统能不能上线的下限很多人把精力都花在研究算法上结果第一步光照不均就把所有后续步骤坑了。我自己的经验是预处理做好后面用最简单的方法都能得到不错的结果预处理没做好再高级的分割算法也救不回来。2.1 灰度化与背景不均校正显微照片最常见的两类输入模式一类是暗颗粒亮背景另一类是亮颗粒暗背景比如部分电镜图。处理的第一步先是统一转成灰度图img imread(particles_sample.jpg); gray rgb2gray(img); % RGB 转灰度已经是灰度的图不影响如果彩色图里颗粒和背景靠颜色区分而亮度接近直接用 rgb2gray 可能丢失对比度。这种情况可以挑对比最强的那个颜色通道做处理比如gray img(:,:,2)或者转成 Lab 色彩空间后用 a、b 通道做阈值分割。灰度图拿到手先别急着分割先用imshow和imhist看一眼灰度分布。我处理过的样品里至少有三分之一存在明显的背景亮度不均中间亮、四周暗或者一侧有阴影。这种情况下直接设一个全局阈值角落里的小颗粒会整体“隐身”。背景不均的校正常规做法是估计背景面然后从原图里减掉它bg imopen(gray, strel(disk, 60)); % 半径要明显大于最大颗粒 corrected gray - bg;imopen是形态学开运算它会在保留小于结构元素的亮细节的同时把大范围的背景亮度变化估计出来。这个结构元素半径的选取是关键必须大于样品里最大颗粒的半径否则开运算会把大颗粒本身也当成背景的一部分。如果颗粒尺寸跨度特别大用imgaussfilt(gray, 50)做高斯模糊估计背景面效果通常也不错。2.2 阈值分割Otsu 只是起点灰度校正之后大多数情况下直方图会呈现双峰形态一个峰是背景一个峰是颗粒。这时候用 Otsu 全局阈值就很合适level graythresh(corrected); bw corrected level; % 暗颗粒亮背景注意反转的问题。graythresh返回的 level 是把图像分成“暗”和“亮”两类的最优阈值至于颗粒到底属于哪一类取决于你的成像条件。暗颗粒亮背景用corrected level亮颗粒暗背景用corrected level。如果直方图不是双峰而是多峰甚至单峰Otsu 的效果会明显变差。多峰往往意味着图像里除了颗粒和背景还有第三种灰度对象比如染色剂残留、抛光布碎屑、反光点。这时候我倾向于在直方图谷底手动定阈值。用imhist拉出灰度分布找到谷底对应的灰度值直接写在代码里bw corrected 85; % 手动阈值85 根据直方图谷底确定还有一种情况是背景校正后仍然存在局部亮度波动比如颗粒本身反光强烈那么固定阈值会在颗粒内部留下空洞。这时用自适应阈值比全局阈值更稳妥bw imbinarize(corrected, adaptive, Sensitivity, 0.4, ForegroundPolarity, dark);Sensitivity 参数从 0 到 1越大越容易把暗区判定为前景。我建议刚开始用默认的 0.4然后结合视觉效果慢慢调。2.3 形态学清理与边缘颗粒取舍分割出来的二值图绝不会是干干净净的。细小噪点、颗粒内部的孔洞、颗粒边缘的毛刺都会干扰后续统计所以要按固定顺序做一轮清理bw bwareaopen(bw, 50); % 删除小于 50 像素的噪点 bw imclose(bw, strel(disk, 3)); % 闭合细小缝隙 bw imfill(bw, holes); % 填充颗粒内部孔洞 bw imclearborder(bw); % 删除触碰到图像边界的连通域这里每一步都有一个值得想清楚的“为什么”。bwareaopen的最小面积阈值要根据你样品里的最小真实颗粒来定过滤的只是比真实颗粒小一个量级的噪点否则会把真颗粒误删。imclose的结构元素半径宁可小一点也不要大因为每膨胀一圈颗粒周长就会被系统性放大。imfill填充孔洞要谨慎如果颗粒本身是疏松多孔的比如某些催化剂载体内部孔洞是有物理意义的填不填取决于你要统计的是“颗粒实体面积”还是“颗粒外轮廓面积”。最稳妥的做法是同一批样品用同一个规则并在报告里注明。imclearborder这一步要单独说一下。颗粒只要有一部分在图像边缘外它的面积、周长、粒径就全是残缺数据。如果你想保证统计精度就删掉它们如果这批样品的大颗粒本来就容易压在边缘删掉后会把粒度分布里的粗端截掉那就要在报告里说明这个偏差。我通常的做法是默认删除因为残缺颗粒对粒度分布是明显的负向污染比少统计几个颗粒更危险。3. 粘连颗粒的拆分与计数分水岭的正确打开方式预处理干净的二值图颗粒之间往往是挨着的。尤其粉末样品压片、制样过程中颗粒天然会堆叠拍出来就是一片粘连的团块。这种粘连不处理计数和测量都会出大问题。3.1 连通域计数粘连会低估数量最基础的计数是连通域分析cc bwconncomp(bw, 8); % 8 邻域连通 numObjects cc.NumObjects; % 颗粒数代码本身没问题问题在于语义bwconncomp把相互接触的前景像素归为同一个连通域。两颗颗粒只要碰在一起就被算成一个大连通域数量被合并面积被叠加周长也被算成了两个颗粒的“合体轮廓”。所以连通域计数只适用于颗粒完全分离的理想样品。现实里粘连比例超过百分之十的样品很常见。这时候就必须在计数和测量之前先把粘连颗粒切开。3.2 标记控制的分水岭拆分流程分水岭算法是处理粘连颗粒最经典的方案。它的直观理解是把二值图里的每个颗粒当作一座山距离背景越远的点“海拔”越高然后模拟水从山顶往山脊漫延不同山头的积水在交界处相遇就形成了分割线。直接对二值图跑分水岭会产生严重的过分割因为颗粒表面的毛刺和噪声都会形成假的山谷。所以正确的做法是标记控制的分水岭流程如下D bwdist(~bw); % 距离变换每个前景像素到最近背景的距离 D -D; % 取负把颗粒内部的中心变成“盆地” mask imextendedmin(D, 2); % 找到深度不小于 2 像素的局部极小值作为种子 D2 imimposemin(D, mask); % 强制只在种子位置产生极小值 L watershed(D2); % 分水岭分割 bw(L 0) 0; % 将分水岭脊线置为背景实现颗粒拆分imextendedmin的第二个参数 H 是关键中的关键。它表示只保留深度不小于 H 像素的局部极小值作为种子。H 太小颗粒表面的微小起伏都会成为独立种子导致一个颗粒被切好几瓣H 太大真实颗粒的种子被过滤掉粘连又分不开。针对大部分粒径在 20 到 200 像素的样品H 从 1 到 5 之间试观察拆分效果的稳定区间。拆分之后重新跑一遍bwconncomp和后续测量得到的数量和每个颗粒的独立面积、周长才是可信的。3.3 两个最常踩的坑过分割与无法拆分第一个坑是过分割。细长颗粒特别容易被分水岭拦腰切断。原因是细长颗粒的距离变换会在长轴方向形成一条“脊”这条脊上的局部极小值点可能不止一个。解决方法是适当调大 H或者对分割后的连通域做一次形状判断把长径比异常高的小碎块合并回去。第二个坑是无法拆分。如果两颗颗粒是大面积平面接触而不是像两个圆球那样只有一个切点它们之间就不会形成明显的距离变换“谷”分水岭没有依据去切分再调参数也没用。这类高粘连样品我的处理策略是先尝试分水岭再看拆分后的统计结果是否合理如果分不开就放弃逐颗粒测量改用下面的面积估算法——假设从样品里随机挑出一批能分离的颗粒统计出它们的平均等效面积 A_mean然后总颗粒数 ≈ 总前景像素面积 / A_mean。这是个粗糙但必要的兜底方案适合粘连比例超过一半、形态学手段几乎失效的情况。4. 周长、面积和粒径的量测regionprops 背后的算法逻辑测量部分我统一走regionprops但这个函数内部的实现细节直接决定了结果精度尤其是周长一项非常值得掰开讲清楚。4.1 面积像素求和与误差边界面积是整条链路里最可靠的一项。二值图里每个像素如果不是前景就是背景连通域的“面积”在像素意义上就是像素个数乘以单个像素的面积。regionprops返回的Area就是像素个数stats regionprops(cc, Area, Perimeter, EquivDiameter);它的统计本质是在数像素没有复杂的插值或拟合。误差主要来自边界像素的归属问题——一粒像素被判定为前景还是背景取决于阈值。给定阈值后像素计数是确定性的灰度图像上颗粒边界在亚像素尺度的信息已经丢失这部分误差属于原理性误差只能靠提高分辨率来缩小。对于粒径在 100 像素以上的颗粒面积测量的相对误差通常能控制在 1% 上下对于只几十像素的小颗粒误差就到 5% 甚至更高了。一个容易被忽略的点如果我直接把二值图传给regionprops它会默认把颗粒内部的孔洞也计入连通的区域属性导致面积偏大。加上imfill这一步保证送入regionprops的都是实心连通域。4.2 周长为什么“数边缘像素”会高估周长的算法就没那么直白了。最简单的想法是提取边界像素然后数个数bwperim求和但这会带来显著误差。原因是斜向相邻的两个边界像素真实距离是 √2 像素而直接数个数会把它们当成 1 像素。一条 45 度倾斜的直线边界用数边界像素的方法算出来的周长比真实值低估了约 30%。regionprops的Perimeter字段做了改进它沿边界逐点追踪水平或垂直相邻的两点贡献 1对角相邻的两点贡献 √2累加起来。这已经是相当合理的估计。不过如果你要更高精度的结果可以自己用链码方式算一次边界长度bounds bwboundaries(bw, 8, noholes); perimTotal 0; for k 1:numel(bounds) b bounds{k}; dx diff(b([1:end 1], 1)); dy diff(b([1:end 1], 2)); perimTotal perimTotal sum(sqrt(dx.^2 dy.^2)); end这段代码把闭合边界上相邻像素中心的欧氏距离累加起来比regionprops的折线近似更接近真实轮廓的“像素中心路径”。需要特别提醒的是周长测量的分辨率敏感性。颗粒边界越放大越能看到毛刺般的细节而每多出一个锯齿周长就增加一段。这跟真实世界里的海岸线测量悖论本质相同。所以跨批次的周长比较必须保证图像分辨率、放大倍数、分割参数完全一致否则测出来的“周长差异”可能只是分辨率差异。4.3 用 regionprops 一次取齐全部指标我实际用的regionprops调用会一次取回全部需要的指标并且同时传入灰度图以便顺便统计每个颗粒的平均灰度stats regionprops(cc, gray, Area, Perimeter, EquivDiameter, ... MajorAxisLength, MinorAxisLength, Orientation, ... FeretDiameter, Centroid, BoundingBox, MeanIntensity);各字段与物理含义对照如下字段含义注意事项Area面积像素数需乘 scale² 换为物理面积Perimeter周长像素路径长度对角步按 √2 计分辨率影响大EquivDiameter等效圆直径√(4A/π)从面积换算而来MajorAxisLength / MinorAxisLength等效椭圆长轴、短轴基于二阶矩旋转不敏感Orientation等效椭圆主轴角度单位是度FeretDiameter最大费雷特径依赖旋转朝向Centroid形心坐标可用于定位统计BoundingBox外接矩形判断是否切边MeanIntensity颗粒区域的平均灰度需要传入灰度图拿到这些字段后圆形度需要自己算C 4*pi*Area ./ (Perimeter.^2)。我习惯把它作为质量指标输出如果一批样品的圆形度普遍在 0.8 以上那说明颗粒接近球形如果普遍低于 0.4就得考虑是不是分割环节出了问题比如过分割或欠分割把单个颗粒扭曲了。5. 像素标定没有这一步所有测量值都是“像素数”在显微镜下图像里的长度没有物理意义完全取决于放大倍数。不标定的话你得到的所有尺寸都停留在“像素”单位没法跟筛分粒度、鲍氏比表面积这些工程数据对上。标定是整个系统里最简单但又最容易被忽略的环节。5.1 三种常用标定方法第一种是标定板法也是精度最高的做法。用与样品相同的显微镜和放大倍数拍一张标定片标定片上刻有已知间距的刻度线。读出刻度线之间的像素距离就能算出每个像素对应的物理长度。第二种是参考物法。在样品里放一个已知直径的物体比如标准微球、金属丝、硬币拍完图之后量出参考物的像素直径。这个方法精度低于标定板但胜在方便适合无标定板的场合。第三种是读图面信息法。很多电镜图自带比例尺或者软件采集界面上会显示视野宽度。比如某个放大倍数下视野宽度是 1000 微米图像宽度是 2048 像素那比例尺就是 1000/2048 微米每像素。这个方法精度取决于厂商标称的准确性用于筛查足够用于精确对比就不够了。无论用哪种方法处理逻辑是一样的在一张参考图像上量出已知长度的像素数算出比例尺pixelLen 150; % 参考长度对应的像素数例如标定刻度间距离 physicalLen 0.1; % 参考长度的物理值单位 mm scale physicalLen / pixelLen; % 每像素对应的 mm 数5.2 换算公式与误差传递拿到 scale 之后所有指标统一换算area_mm2 stats.Area .* scale.^2; perimeter_mm stats.Perimeter .* scale; diameter_mm stats.EquivDiameter .* scale;这里有个误差放大关系需要心里有数长度量纲的误差是一阶放大面积量纲的误差是二阶放大。也就是说假如标定板读数有 1% 的误差那么粒径指标有约 1% 的系统误差面积指标却会有约 2% 的系统误差。所以标定这一步值得重复读数取平均而不是随意量一次。我一般会在标定图上选择间隔不同距离的刻度测 3 到 5 组取 scale 的中位数。5.3 一套必须坚持的标定纪律我踩过的最大一次坑是换了物镜之后忘了重新标定拿旧的比例尺去套新图像导致整批粒径结果整体偏了百分之四十多。这个错误的隐蔽性在于数据看起来依然“合理”只是全部偏大或偏小如果不做人工校验根本发现不了。从那以后我定了三条纪律第一每次改变物镜或变倍倍数必须重新标定并做记录第二比例尺和每张图的元数据文件名、物镜倍数、标定日期写在同一张结果表里第三同批次对比一律使用同一套比例尺。宁可数据少不可标定乱。6. 完整系统代码从单图处理到批量报表把前面所有环节拼起来就是一套可以从单图扩展到批量目录的完整系统。我习惯把处理核心封装成一个函数输入图像路径和参数输出一张统计表。6.1 核心函数封装function T analyzeParticles(imgPath, scale, params) % analyzeParticles 对单张图像执行颗粒测量 % scale : 每像素对应的物理长度(mm) % params : 结构体包含 minArea, H, thresholdMode 等 img imread(imgPath); if size(img,3) 3 gray rgb2gray(img); else gray img; end % 1. 去噪与背景校正 gray medfilt2(gray, [3 3]); bg imopen(gray, strel(disk, params.bgRadius)); corrected gray - bg; % 2. 阈值分割 if strcmp(params.thresholdMode, otsu) level graythresh(corrected); bw corrected level; else bw corrected params.customLevel; end % 3. 形态学清理 bw bwareaopen(bw, params.minArea); bw imclose(bw, strel(disk, 2)); bw imfill(bw, holes); bw imclearborder(bw); % 4. 可选分水岭拆分 if params.doWatershed D bwdist(~bw); D -D; mask imextendedmin(D, params.H); D2 imimposemin(D, mask); L watershed(D2); bw(L 0) 0; end % 5. 连通域与指标测量 cc bwconncomp(bw, 8); stats regionprops(cc, gray, Area, Perimeter, EquivDiameter, ... MajorAxisLength, MinorAxisLength, Orientation, ... FeretDiameter, Centroid, BoundingBox, MeanIntensity); T struct2table(stats); % 6. 物理单位换算与派生指标 T.Area_mm2 T.Area * scale^2; T.Perimeter_mm T.Perimeter * scale; T.EquivDiameter_mm T.EquivDiameter * scale; T.Circularity 4*pi*T.Area ./ (T.Perimeter.^2); end这个函数把参数全部收敛到params结构体里好处是调参时不需要改动主流程代码只需要更新参数。我也强烈建议把分割结果另存为一张检查图每处理一张图都能看到“它到底把颗粒抠成什么样了”否则一个隐藏的错误参数可能污染整批结果。6.2 批量文件夹处理与结果合并单图不是目的真实需求几乎都是“这一批 50 张图给我出个汇总”。批量循环的写法很简单files dir(samples/*.jpg); allData []; for i 1:numel(files) imgPath fullfile(files(i).folder, files(i).name); T analyzeParticles(imgPath, scale, params); T.ImageName repmat({files(i).name}, height(T), 1); allData [allData; T]; % 累计合并所有统计行 end这里有个性能提示单张图的耗时主要花在imread和regionprops上。如果你处理的是几千张高分辨率图像逐张循环会比较慢可以改用parfor并行循环同时对输出文件的表头做一次一致性拼接。我个人不建议为了提速去降低图像分辨率因为粒度分布里的小颗粒端对分辨率最敏感。6.3 粒径分布图与 Excel 输出汇总数据之后最常用的输出是粒径分布直方图figure; histogram(allData.EquivDiameter_mm, 25); xlabel(等效圆直径 (mm)); ylabel(颗粒数量); title(粒径分布直方图);直方图的分组数会影响分布形状25 个柱是我常用的起点。更严谨的做法是用固定柱宽这样不同批次的图可以直接叠图对比而不是因为分组数不同导致“峰型变了”的错觉。数据落地用writetable直接写 Excel方便下游人员继续做统计writetable(allData, particle_stats.xlsx);字段里既有像素原始值又有物理单位值这样后续如果发现标定错误还可以用原始像素值重新换算不用重新跑一遍图像。7. 实测踩坑记录与参数调优心法系统能跑通是一回事能稳定输出正确结果是另一回事。以下这几个问题基本是我在实际样品上遇到的问题里复现率最高的。7.1 光照不均让小颗粒“隐身”有一次处理一批白色粉末的显微照片背景是浅灰色的载物台颗粒是半透明的。全局阈值下视野中央的大颗粒全部正常分割但角落里的细颗粒几乎全部丢失。最初我以为是阈值没设好试来试去发现是光照不均角落亮度偏低细颗粒的灰度已经跟背景混在一起了。解决方法是先做背景校正再做阈值而且这一步必须在分割前完成不能偷懒。具体实现就是第 2 节的imopen背景估计。另外medfilt2去噪一定要放在背景校正之前否则中值滤波会把颗粒边缘的微小灰度变化抹平影响后续分割。7.2 颗粒切边过滤还是保留切边颗粒是个两难问题。过滤掉会损失大颗粒样本让分布曲线左偏保留下来残缺颗粒的面积、周长、粒径全是错的会让分布曲线左右都乱。我的处理策略是默认imclearborder全部过滤然后在报告里单独注明“已剔除边界颗粒”。这样统计口径清楚批间可对比。如果你担心大颗粒被过度剔除那就不做边界剔除而是把BoundingBox靠近图像边缘的记录标注出来在汇总时给出一个剔除边界记录的数量和比例由下游决策者自行选择。7.3 高粘连样品分水岭解决不了怎么办有一类样品颗粒在制样时就被压成致密的团块彼此之间根本没有清晰边界。分水岭在这种情况下不仅无效有时候还会因为生成大量假分割线把连通域切得支离破碎。对这种样品我的经验是分清主次如果颗粒形态和材料本身已经确定只是想知道粒度范围那就用第 3.3 节说的面积估算法重点做全视野总面积的统计。如果必须逐颗粒测量那问题已经不属于单纯图像处理的范畴需要回到制样环节解决——分散得更开、或改用动态图像法这类专门的粒度仪。图像处理不是万能的承认某个环节的局限反而能让整个系统更可靠。7.4 验证系统准确性拿人工计数当裁判系统写完之后必须做一次验证才能投入使用。方法很朴素随机挑 5 个视野人工在图像上逐一标记计数再和系统输出对比。误差在 5% 以内说明分割逻辑基本可靠如果系统数得明显偏多多半是过分割数得偏少多半是粘连没有拆开。这个验证不只是做一次换新样品类型、新成像条件时都要重做。我通常还会做一次“参数敏感性测试”把阈值、H、minArea 这三个参数分别往两个方向调 20%看统计结果波动多大。如果波动在可接受范围内说明参数选择落在稳定区间系统对参数不敏感可以放心批量跑如果某一个参数稍微一动结果就大幅变化说明当前工作点太脆弱需要重新找稳定区间。整套系统做到这一步已经不只是个“测粒径的脚本”而是一套能反复验证、标定可控、带质量审计的测量工具。回头再看最花时间的从来不是写代码而是搞清楚每种参数改变之后分割结果到底在图片层面发生了什么变化。所以每处理一批新样品我都会把中间产生的二值分割图翻出来看一遍这比盯着小数点后四位有没有跳更值得。