图像分割这个方向大家听得最多的可能是深度学习那一套U-Net、Mask R-CNN、Transformer 分割头动辄要标注几千张图、训上几天几夜。但实际工程里尤其是医学影像、工业检测、材料显微分析这些场景你根本拿不到标注数据或者标注成本高到离谱。这时候无监督方法就成了刚需。局部质心分割就是其中一类思路非常朴素、但实战里意外好用的方法——它不需要训练、不需要标注靠像素或体素在空间上的密度分布自动找出图像里的团块中心再据此把前景和背景掰开。2D 图像能用3D 体数据同样能用而且核心代码用 Matlab 写出来也就百来行。这篇就把这套方法的原理、Matlab 实现、参数调优和踩坑经验完整拆一遍适合做医学图像、点云切片、工业缺陷检测以及任何没标注但想先跑出个分割结果看看的朋友。1. 局部质心分割到底在解决什么问题1.1 从找中心这个直觉说起先别急着看公式。你想象一张细胞显微图背景是均匀的浅灰中间散布着若干个深色的细胞团。人眼怎么区分前景背景其实靠的是哪里密集、哪里稀疏。细胞团内部像素扎堆背景像素稀稀拉拉。局部质心方法就是把这个直觉数学化在图像上滑动一个小窗口计算窗口内像素的加权重心这个重心就是局部质心。如果某个像素离它所在窗口的质心很近说明它处在一个密集区域大概率是前景如果离得很远说明它处在稀疏区域大概率是背景。这个思路和 K-means、Mean Shift 有亲缘关系但比它们轻量得多。K-means 要预设簇数量Mean Shift 要迭代收敛而局部质心法基本是一遍扫描加一个阈值判断计算量小、无需迭代、无需预设类别数。这就是它在无监督场景下的最大优势——你不需要知道图里有几个目标也不需要告诉算法这是几分类。1.2 无监督这三个字在工程里的分量很多人做项目时对无监督有误解以为无监督就是效果差一点但省事。实际恰恰相反在标注数据稀缺的领域无监督方法往往是唯一能落地的方法。比如医学图像分割一个资深医生标注一张 3D CT 要花几十分钟一个项目动辄几百例标注预算直接爆炸。再比如工业流水线上的缺陷检测缺陷形态千变万化你标注了划痕明天来了个气泡模型就懵了。无监督方法不依赖具体形态只依赖前景比背景更聚集这个统计假设泛化性反而更好。局部质心法的适用前提很明确前景目标在空间上是成团的、密度高于背景。只要满足这个前提它就能工作。不满足的场景比如前景是细长的线状结构、或者前景背景密度差不多那它就会失效。这一点必须在动手前想清楚否则后面调参调到怀疑人生。1.3 2D 和 3D 为什么能用同一套逻辑2D 图像里窗口是个矩形质心是 (x, y) 两个坐标。3D 体数据里窗口是个立方体质心是 (x, y, z) 三个坐标。数学形式完全一样只是维度加了一维。这意味着同一套 Matlab 代码稍微改一下窗口定义和坐标计算就能从 2D 直接迁移到 3D。对于做医学影像CT/MRI 是 3D、材料断层扫描也是 3D的朋友这个迁移成本几乎为零非常划算。提示3D 数据的计算量是 2D 的立方级增长。一个 512×512 的 2D 图窗口 15×15 也就 225 个像素换成 512×512×200 的体数据窗口 15×15×15 就是 3375 个体素再乘以体素总数内存和耗时都要提前评估。2. 局部质心的数学原理与窗口设计2.1 质心计算的完整推导设图像为 $I$在位置 $p$ 处取一个半径为 $r$ 的邻域窗口 $W_p$。窗口内每个像素有一个权重最简单的权重就是像素灰度值本身灰度越高权重越大也可以先用边缘检测或梯度幅值做权重。局部质心 $c_p$ 定义为$$c_p \frac{\sum_{q \in W_p} w(q) \cdot q}{\sum_{q \in W_p} w(q)}$$其中 $q$ 是窗口内像素的坐标$w(q)$ 是权重。这个公式说白了就是加权平均位置。如果窗口内所有像素权重一样质心就是窗口的几何中心如果窗口内有一团高权重像素偏在一边质心就会被拉过去。得到质心后计算当前像素 $p$ 到质心的距离 $d_p |p - c_p|$。这个距离就是局部质心距离特征。前景像素因为周围都是同类高权重像素质心基本就在自己附近$d_p$ 小背景像素周围稀疏质心被拉向远处的前景团$d_p$ 大。于是一个简单的阈值就能把两者分开。2.2 窗口半径 r 的选择逻辑窗口半径是整个方法里最关键的参数没有之一。$r$ 太小窗口内像素太少质心估计噪声大分割结果碎成一片$r$ 太大窗口跨越了前景和背景质心被平均掉前景背景区分度下降。经验上$r$ 应该略大于目标特征的半宽度。比如细胞直径约 30 像素那 $r$ 取 15 到 20 比较合适。实际操作中我一般这样定先目测或统计前景目标的典型尺寸取 $r$ 为目标半径的 1.0 到 1.5 倍然后在这个范围里做小步长扫描比如 10、12、15、18、20看哪个 $r$ 下质心距离的直方图双峰最明显。双峰越明显说明前景背景可分性越好。这个方法比拍脑袋靠谱得多。2.3 权重函数的选择与影响权重函数决定了什么算密集。最朴素的是用原始灰度但灰度受光照不均影响大。更稳的做法有两种一是用梯度幅值的倒数边缘处权重低平坦区域权重高二是先做局部对比度归一化再取灰度。我在工业检测项目里更倾向第二种因为金属表面反光会导致灰度剧烈波动归一化之后质心估计稳定很多。还有一种进阶玩法用局部方差做权重。前景内部方差小纹理一致背景方差大噪声杂乱用方差的倒数当权重能进一步拉开前景背景的差距。这个技巧在低信噪比的荧光显微图里效果特别明显。3. Matlab 实现从 2D 到 3D 的完整代码拆解3.1 2D 版本的核心函数先看 2D 的实现。核心就是两层循环加一个窗口求和但直接写循环在 Matlab 里慢得让人抓狂必须向量化。我的做法是用imfilter或者conv2一次性算出所有位置的加权坐标和。function [seg, distMap] localCentroid2D(I, r, wMode) % I: 输入灰度图, double 类型, 归一化到 [0,1] % r: 窗口半径 % wMode: 权重模式 gray | grad | var if nargin 3, wMode gray; end I im2double(I); [H, W] size(I); % 构造权重图 switch wMode case gray Wt I; case grad [gx, gy] gradient(I); Wt 1 ./ (sqrt(gx.^2 gy.^2) 1e-3); case var mu imfilter(I, fspecial(average, 2*r1), replicate); mu2 imfilter(I.^2, fspecial(average, 2*r1), replicate); v max(mu2 - mu.^2, 1e-6); Wt 1 ./ v; end % 坐标网格 [X, Y] meshgrid(1:W, 1:H); % 加权坐标和用均值滤波实现窗口求和 k fspecial(average, 2*r1); sumW imfilter(Wt, k, replicate); sumWX imfilter(Wt .* X, k, replicate); sumWY imfilter(Wt .* Y, k, replicate); % 局部质心 Cx sumWX ./ (sumW 1e-8); Cy sumWY ./ (sumW 1e-8); % 质心距离图 distMap sqrt((X - Cx).^2 (Y - Cy).^2); % Otsu 自动阈值 lvl graythresh(mat2gray(distMap)); seg distMap lvl * max(distMap(:)); end这段代码的关键点在于用均值滤波代替显式窗口循环。imfilter配合fspecial(average, 2*r1)本质上就是窗口求和Matlab 底层是优化过的比手写 for 循环快几十倍。sumWX和sumWY分别算的是窗口内加权 x 坐标和 y 坐标之和除以权重和就得到质心坐标。整个函数没有一处显式循环跑一张 1024×1024 的图也就几十毫秒。3.2 3D 版本的迁移改动3D 版本改动很小主要是把fspecial换成fspecial3或者自己构造 3D 均值核坐标网格加一维。function [seg, distMap] localCentroid3D(V, r, wMode) V im2double(V); [H, W, D] size(V); switch wMode case gray Wt V; case grad [gx, gy, gz] gradient(V); Wt 1 ./ (sqrt(gx.^2 gy.^2 gz.^2) 1e-3); case var k ones(2*r1,2*r1,2*r1) / (2*r1)^3; mu imfilter(V, k, replicate); mu2 imfilter(V.^2, k, replicate); Wt 1 ./ max(mu2 - mu.^2, 1e-6); end [X, Y, Z] ndgrid(1:H, 1:W, 1:D); k ones(2*r1,2*r1,2*r1) / (2*r1)^3; sumW imfilter(Wt, k, replicate); sumWX imfilter(Wt .* X, k, replicate); sumWY imfilter(Wt .* Y, k, replicate); sumWZ imfilter(Wt .* Z, k, replicate); Cx sumWX ./ (sumW 1e-8); Cy sumWY ./ (sumW 1e-8); Cz sumWZ ./ (sumW 1e-8); distMap sqrt((X-Cx).^2 (Y-Cy).^2 (Z-Cz).^2); lvl graythresh(mat2gray(distMap)); seg distMap lvl * max(distMap(:)); end注意 3D 里ndgrid和meshgrid的区别ndgrid生成的维度顺序和数组本身一致用imfilter时不会错位。我早期用meshgrid踩过坑前两维对得上第三维坐标全乱了分割结果在 z 方向完全错位排查了半天才发现是网格函数用错了。3.3 内存与速度的实测数据拿一组实测数据说话。测试环境是 16GB 内存的普通工作站Matlab R2023b。数据类型尺寸窗口半径耗时峰值内存2D 灰度1024×1024150.08 s约 120 MB2D 灰度2048×2048150.35 s约 480 MB3D 体数据256×256×12872.1 s约 1.8 GB3D 体数据512×512×200718 s约 9 GB可以看到 3D 的内存增长非常凶。512×512×200 的体数据光distMap一个 double 数组就是 512×512×200×8 字节 ≈ 420MB再加上sumWX、sumWY、sumWZ、Wt这些中间变量轻松上 9GB。所以 3D 场景下我强烈建议先把体数据降采样到 256 级别跑通流程确认参数后再上原始分辨率或者分块处理。4. 参数调优与分割质量提升的实战技巧4.1 阈值选取Otsu 不是万能的上面代码里用了graythreshOtsu 法自动定阈值。Otsu 在质心距离直方图呈明显双峰时很好用但实际数据往往不是理想双峰。比如前景目标大小差异很大小目标的质心距离分布和大目标重叠Otsu 就会切错。我的经验做法是先画质心距离的直方图看一眼。如果双峰明显Otsu 直接用如果是一个主峰带一个长尾那就用百分位数阈值比如取距离图的 70% 分位数作为阈值。这个 70% 不是固定的要根据前景占比预估。如果你大概知道前景占 20%那阈值就取 80% 分位数附近。这个先看直方图再定阈值的习惯帮我省下了大量盲目调参的时间。4.2 后处理形态学开闭运算的取舍原始分割结果几乎一定会有噪点和小孔洞。标准做法是开运算去噪点、闭运算填孔洞。但这里有个坑开运算的核大小不能超过目标的最小尺寸否则小目标直接被抹掉。我一般用半径 2 到 3 的圆盘核做开运算闭运算核稍大一点半径 3 到 5。对于 3D 数据形态学操作要用 3D 结构元素。Matlab 里strel(sphere, r)生成球形结构元素比立方体结构元素更符合各向同性的物理直觉。但要注意如果体数据的层间距z 方向分辨率和层内分辨率差很多球形核会被拉伸这时候应该用椭球核z 方向半径按层间距比例缩放。4.3 多尺度融合解决目标尺寸差异单一窗口半径只能适配一种目标尺寸。如果图里既有大团块又有小颗粒单尺度必然顾此失彼。解决办法是多尺度用 3 个不同半径比如 r、2r、4r分别算质心距离图然后取每个像素在三个尺度下的最小值作为最终特征。因为不管目标多大总有一个尺度能匹配上它那个尺度下的质心距离就会很小。这个多尺度融合的思路在细胞分割里特别有效因为细胞核大小差异经常有两三倍。代价是计算量翻三倍但换来的是分割完整度的大幅提升值。注意多尺度融合后阈值选取要重新做。因为最小值操作会把整体距离分布压低原来的阈值不再适用必须重新看直方图。5. 踩坑实录那些让我熬夜的失败案例5.1 光照不均导致的质心漂移早期做一个金属表面缺陷检测项目图像一侧有强反光另一侧偏暗。直接用灰度当权重亮的那侧质心全被拉到反光区暗侧缺陷完全分不出来。当时我以为是窗口半径没调好试了十几个 r 值都没用后来才意识到是权重的问题。解决办法是先用imflatfield或者大核背景减除去掉低频光照分量再做质心计算。Matlab 的imflatfield在 R2018b 之后就有一行代码搞定。这个坑的教训是局部质心法对光照不均极其敏感因为光照直接改变了权重分布。任何无监督方法预处理的重要性都不亚于算法本身。5.2 3D 数据的各向异性陷阱医学 CT 数据经常是层内 0.5mm 分辨率、层间 2mm 分辨率各向异性比达到 4:1。我一开始直接用球形窗口结果 z 方向的分割边界糊成一片。原因是球形窗口在物理空间里其实是个被压扁的椭球z 方向实际覆盖的物理距离是层内的 4 倍质心估计在 z 方向严重平滑。修正方法是把窗口定义在物理坐标系里而不是体素坐标系。具体做法是给 z 方向的坐标乘以层间距比例或者直接用椭球窗口z 半径设为 r 除以各向异性比。这个细节在大部分教程里都不会提但做医学影像的人一定会遇到。5.3 边界目标的截断效应图像边缘的目标窗口会超出图像范围。如果直接补零边缘处权重和骤降质心估计完全失真边缘目标会被错误分割。正确做法是用replicate边界模式代码里已经用了让窗口超出部分复制边缘像素。但即便如此图像最外圈几个像素的分割结果仍然不可靠实际使用时要裁掉边缘一圈。5.4 计算精度double 还是 singleMatlab 默认 double但 3D 大数据用 double 内存吃不消。我试过改 single结果质心距离在阈值附近出现抖动分割边界变得毛糙。原因是 single 只有 7 位有效数字而质心坐标计算涉及大数相减坐标值减去质心值有效数字损失严重。折中方案是权重图和坐标和用 double 累加最终距离图转 single 存储。这样内存省一半精度基本不受影响。6. 典型应用场景与效果对比6.1 医学影像CT 肺结节初筛肺结节在 CT 里是高密度小团块背景是低密度肺实质密度差异明显非常适合局部质心法。我用 3D 版本在公开肺结节数据上跑过窗口半径按结节典型半径 5mm 设定配合多尺度融合召回率能到 85% 以上。虽然精度不如训练好的深度学习模型但胜在零标注、零训练作为初筛工具完全够用能把医生需要看的候选区域从整个肺缩减到几十个。6.2 工业检测显微图像颗粒计数材料显微图里的颗粒、粉末、夹杂物都是成团的高密度区域。局部质心法不仅能分割还能顺便做计数——分割后的连通域个数就是颗粒数。这里的关键是窗口半径要匹配最小颗粒尺寸否则小颗粒会被漏掉。我一般取最小颗粒半径的 1.2 倍作为 r宁可稍微过分割也不要漏检因为后续可以用面积阈值过滤掉过分割的碎片。6.3 与深度学习方法的定位差异必须说清楚局部质心法不是要取代 U-Net 这类方法它们的定位完全不同。深度学习需要标注、需要训练、需要 GPU但精度高、能处理复杂语义。局部质心法零成本、可解释、能直接跑在 CPU 上但只能处理密度可分的简单场景。实际项目里我经常把两者结合先用局部质心法快速生成候选区域和伪标签人工修正一小部分后拿去训练深度学习模型。这个无监督预筛 弱监督精修的流程能把标注成本降低一个数量级。对比维度局部质心法深度学习分割标注需求无大量训练时间无数小时到数天硬件需求CPU 即可通常需 GPU可解释性强弱复杂语义处理不支持支持适用场景密度可分团块通用7. 代码工程化与批量处理的几点建议7.1 封装成可配置的函数上面给的函数是最小可用版本。实际项目里我会把它封装成一个类把窗口半径、权重模式、阈值方法、后处理参数都做成属性这样调参时不用改代码改配置就行。Matlab 的arguments块R2019b 之后做参数校验很方便能避免传错类型导致的诡异错误。7.2 批量处理与并行加速处理几十上百例数据时用parfor并行是最直接的加速手段。但要注意parfor里不能有跨迭代的共享状态每个 worker 独立处理一例数据。另外 3D 数据内存占用大并行数不能开太多否则内存爆掉。我的经验是并行数设为min(核心数, floor(可用内存 / 单例峰值内存))宁可少开几个也别爆内存。7.3 结果可视化与质控分割结果一定要可视化检查尤其是 3D 数据。Matlab 的volshowR2019b 之后能直接渲染 3D 分割结果配合labeloverlay在 2D 切片上叠加轮廓质控效率很高。我习惯把原始图、质心距离热图、分割结果三联图输出一眼就能看出问题出在哪个环节——是权重不对、还是阈值不对、还是后处理过度。7.4 参数持久化与复现调好一组参数后一定要连同数据版本、Matlab 版本一起记录下来。我吃过亏同一个项目隔了半年重新跑Matlab 升级了graythresh的默认行为有细微变化结果对不上。现在我的习惯是把所有参数存成 JSON和结果一起归档保证任何时候都能复现。8. 从局部质心延伸出去的几个改进方向局部质心法本身很简单但它的框架可以往上叠很多改进。第一个方向是自适应窗口不再用固定半径而是根据局部密度动态调整窗口大小密集区域用小窗口、稀疏区域用大窗口。第二个方向是迭代精化把第一次分割结果当权重重新算质心迭代两三次边界会明显收敛这其实就是 Mean Shift 的思路。第三个方向是结合边缘信息质心距离负责找团块边缘检测负责定边界两者融合能显著改善边界贴合度。我自己最常用的是迭代精化实现简单、收益明显。一般迭代 3 次就收敛了再多收益递减。代码上就是在原函数外面套一层循环把seg转成权重再喂回去。这个改动不到十行但分割边界的平滑度和准确度都有肉眼可见的提升。最后分享一个我在实际项目里总结的小经验局部质心法的参数没有最优解只有针对当前数据的最优解。换一批数据窗口半径和阈值可能都要重调。所以别指望调好一组参数就一劳永逸把调参流程标准化、把质控可视化做好比追求某个神奇参数值重要得多。这套方法的价值不在于它多精确而在于它能在你什么都没有的时候快速给你一个能看的、能迭代的起点。