k-means彩色图像分割：Matlab源码拆解与多尺度金字塔调优

简介这是一份基于k-means聚类算法的Matlab彩色图像分割项目源码面向图像处理入门及进阶开发者旨在通过聚类实现像素级颜色区域划分适用于目标提取、图像前景/背景分离等常见场景。压缩包共55个文件其中19个.m脚本为核心算法与工具函数涵盖预处理、特征计算、聚类合并、评估等环节另附25张jpg测试图片、若干.mat数据与文档便于直接复现与对比效果整体体积仅869KB。资源在发布前经过亲测校正可稳定运行若遇问题还可联系作者获得指导。已有1207人学习下载口碑较好。整套代码结构清晰、注释完整既能作为课程设计或毕设的参考实现也适合用于理解k-means在图像分割中的实际调参与优化思路。1. 一张彩色照片k-means为何仍是分割的及格线手头有张starfish.jpg海星和海水颜色分明但边缘有一圈渐变光晕。用深度学习做分割要标注、要训练周期按天算用 matlab 自带的kmeans直接丢进去又常常把高光区域单独分出来。这套kmeansClusters项目的价值在于它把彩色图像分割拆成了预处理、距离计算、多尺度聚类、簇合并、质量评估五个环节并且每个环节都有独立函数改一个参数就能看到分割边界变化。对于课程设计、算法新人和需要快速验证特征思路的工程师这套源码比直接调kmeans更有教学意义也比上深度学习模型更容易排查问题。下面从文件结构开始逐层拆解这条分割链路。2. 源码结构拆解从彩色原图到可复现的分割标签2.1 解压之后的文件分类先分清哪些该改、哪些只该读打开资源包第一眼是几十个.m文件和一堆.jpg测试图。如果直接挨个点开很容易迷失。按数据流方向可以把文件分成四组如下表分组关键文件职责测试入口test1.m、jefftest.m、angitest.m调用分割流程并可视化结果核心算法preprocess.m、segment.m、kmeansClusters.m特征提取、流程编排、聚类迭代辅助计算sqdist.m、distances.m、createPyramid.m距离矩阵、金字塔生成、特征转换后处理与评估clusterMerge.m、formCoreClusters.m、calculateQ.m、statsCluster.m合并小簇、定量评价、统计输出测试图命名也有讲究。starfish.jpg、beans.jpg这类前景背景对比度高的K 值设 2 或 3 就能看到清晰分割buildings.jpg、angkorwat2.jpg这类包含天空、建筑、植被、阴影的多区域场景至少需要 5 个簇bagofgrain.jpg、patchofgrass.jpg这类纹理密集的如果不开后处理合并分割结果会像马赛克。2.2 预处理把像素变成聚类算法真正需要的特征矩阵k-means 聚类的输入是特征矩阵不是图像本身。preprocess.m的核心任务就是把图像的每个像素映射成一维特征向量。最简单的方式是取 RGB 三个通道但实际效果往往不好因为 RGB 通道之间存在强相关性且对光照敏感。一个常见的做法是转换到 Lab 颜色空间function [feat, height, width] preprocess(img, useLab) % 将彩色图像转为N*3特征矩阵N为像素数 if nargin 2 useLab true; end if useLab cform makecform(srgb2lab); lab applycform(img, cform); % 使用L、a、b三个通道作为特征 feat double([lab(:,:,1) lab(:,:,2) lab(:,:,3)]); else % 直接用RGB feat double([img(:,:,1) img(:,:,2) img(:,:,3)]); end % 记录原始尺寸 [height, width, ~] size(img); % 转成N*3矩阵 feat reshape(feat, height * width, 3); % 归一化到[0,1]区间 minV min(feat, [], 1); maxV max(feat, [], 1); feat bsxfun(minus, feat, minV); feat bsxfun(rdivide, feat, maxV - minV 1e-6); end这段代码里有几个容易被忽略的地方。第一makecform在 matlab R2014b 之后依然可用但新版推荐使用rgb2lab不过为了兼容老代码源码里保留makecform很正常。第二归一化必须按列计算如果用整个矩阵的全局最大值会导致暗部通道的数值被压得过小。第三reshape的顺序是按列填充所以后面把标签 reshape 回图像尺寸时要使用相同的列主序规则否则分割图和原图会错位。如果图像尺寸太大比如 2000×1500逐像素特征矩阵就有 300 万行直接跑 k-means 会非常慢。这时可以先用imresize降采样或者让preprocess.m返回一个降采样后的尺寸配合后面的金字塔完成多尺度计算。2.3 segment.m 管线把拆分好的零件装成一台机器segment.m是整个项目的门面输入一张图像和参数输出分割标签。它内部的时间线大致如下先做预处理然后判断是否启用金字塔接着从低分辨率层开始聚类最后做簇合并。一个典型的调用方式img imread(starfish.jpg); [labels, centers] segment(img, K, 3, PyramidLevels, 2, DoMerge, true);这里的K是簇数量PyramidLevels是金字塔层数DoMerge决定是否在聚类后合并小簇。注意segment.m的实际实现可能使用输入参数解析结构体pairs所以在命令行用名值对传参时一定要确保参数名和源码里完全一致大小写都不能差。如果出现Undefined function segment或Invalid parameter name多半是路径没加对用addpath把整个资源目录加进来。segment.m里最容易踩坑的是它可能在金字塔最低层使用全部特征维度但在高层只使用颜色特征而丢掉空间坐标。这会导致顶层聚类收敛后映射到底层时出现边缘锯齿。这个问题没有标准解法需要你根据测试图像自行决定是否在特征中加入像素坐标。加入坐标的方式是在preprocess.m输出后拼接归一化的x、y代价是聚类时间变长收益是区域连续性显著提升。3. k-means 核心实现kmeansClusters.m 与 sqdist.m 的距离计算3.1 为什么不用 matlab 内置 kmeans 函数很多人会问matlab 的kmeans成熟稳定为什么要自己写kmeansClusters.m原因有三个。第一内置函数对输入矩阵有各种检查每次调用都产生额外开销在百万级像素点上很吃亏。第二内置函数不方便在每次迭代后插入自定义逻辑比如记录中心移动轨迹、强制合并空簇、观察收敛曲线。第三教学场景下自己写一遍能清楚看到初始化策略、停止条件对分割结果的影响。kmeansClusters.m的大循环分为两步分配assignment和更新update。分配步骤计算每个像素到所有聚类中心的距离把像素归到最近的中心更新步骤重新计算每个簇的质心。下面是一段符合项目风格的精简实现function [labels, centers] kmeansClusters(feat, K, maxIter) % 自定义k-means聚类用于彩色图像分割 % feat: N*D矩阵; K: 簇数; maxIter: 最大迭代次数 N size(feat, 1); % 使用k-means初始化降低随机初始化的影响 centers zeros(K, size(feat, 2)); centers(1, :) feat(randi(N), :); for k 2:K % 计算每个样本到已有中心的最近距离 D sqdist(feat, centers(1:k-1, :)); minDist min(D, [], 2); prob minDist / sum(minDist); cdf cumsum(prob); r rand(); centers(k, :) feat(find(cdf r, 1), :); end labels zeros(N, 1); for iter 1:maxIter % 分配步骤 D sqdist(feat, centers); [~, labels] min(D, [], 2); % 更新步骤 newCenters zeros(size(centers)); for k 1:K idx (labels k); nk sum(idx); if nk 0 % 空簇从样本中随机取一点重新初始化 newCenters(k, :) feat(randi(N), :); else newCenters(k, :) mean(feat(idx, :), 1); end end % 如果中心基本不动提前结束 if norm(newCenters - centers, fro) 1e-6 break; end centers newCenters; end end代码里的 k-means 初始化用cdf累积概率让距离已有中心远的像素被选中为中心的概率更大。这是避免局部最优的关键。如果没有这一步在buildings.jpg这种大面积同色天空的图上两个初始中心可能都落在天空区域另一个远处的建筑区域没有中心最终分割结果会把天空拆成两块。3.2 sqdist.m 的展开式与数值稳定性sqdist.m是kmeansClusters.m里被调用最多的函数它计算两个矩阵之间的平方欧氏距离。朴素实现是三层循环但那在图像像素级场景下不可接受。项目里的实现使用展开式function D sqdist(A, B) % 计算A(NxD)和B(MxD)矩阵中每一对向量的平方欧氏距离 % 返回N*M矩阵DD(i,j) ||A(i,:) - B(j,:)||^2 AA sum(A .^ 2, 2); % N*1 BB sum(B .^ 2, 2); % M*1 D bsxfun(plus, AA, BB) - 2 * (A * B); D(D 0) 0; % 清除浮点误差 end这个展开式的数学依据是||a-b||^2 ||a||^2 ||b||^2 - 2ab利用矩阵乘法一次算出所有点对距离复杂度从三重循环降为矩阵乘法。bsxfun(plus, AA, BB)计算每对向量的模长和A*B是内积矩阵。最后一行D(D 0) 0很重要因为大矩阵乘法里会出现-1e-12这类负值虽然不影响min取索引但会影响后面基于距离的目标函数计算。如果使用的是较新版本的 matlabbsxfun可以替换为隐式扩展AA BB但项目为了兼容 R2016a 之前的版本保留了bsxfun写法。你在二次开发时如果不再关心老版本可以替换速度会略有提升。3.3 空簇、局部最优与迭代停止条件的取舍图像分割场景里空簇出现的概率远高于普通数据集。原因是一张图中某些颜色区域可能只占几百个像素初始化中心如果都集中在大区域小区域很容易在迭代中失去所有像素。处理空簇的方式有两种一种是这里的随机重新初始化另一种是把空簇分裂成当前最大簇的一部分。随机重新初始化实现简单但可能引入噪声分裂大簇能保持几何结构但实现复杂。对于课程设计随机初始化足够用。迭代停止条件除了“中心几乎不动”还可以设置“标签变化比例小于阈值”。标签变化比例更能反映分割结果是否稳定因为中心移动很小但标签可能还在跳变。你可以把停止条件改成changed sum(labels ~ oldLabels) / N; if changed 1e-3 break; end实际经验是maxIter设为 30 到 50在这套测试图上基本都能收敛。设成 100 并不会让效果更好只会让程序在异常数据上多空转两倍时间。4. 多尺度金字塔与 K 值参数实验不同测试图的最优配置4.1 createPyramid.m 降采样后的计算优势segment.m里如果不加金字塔k-means 在buildings.jpg这种 1024×768 的图上要迭代 30 次每次计算 78 万像素到 K 个中心的距离总耗时以分钟计。createPyramid.m先把图像逐层缩小在最顶层只处理几万像素然后把顶层的聚类中心作为下一层的初始中心每一层只需迭代很少次数。createPyramid.m的核心就是 matlab 自带的impyramidfunction pyr createPyramid(img, levels) % 生成彩色图像的高斯金字塔 % levels: 金字塔层数pyr{1}为原图pyr{levels}为最小图 pyr cell(levels, 1); pyr{1} img; for l 2:levels pyr{l} impyramid(pyr{l - 1}, reduce); end end这里有一个隐藏的坑impyramid要求输入图像每个维度至少为 2 的整数次幂关系否则输出尺寸会向下取整导致边缘像素被丢弃。比如一个 1023×700 的图像reduce 一次变成 511×350再 reduce 一次变成 255×175到第 4 层就成了 127×87边缘信息丢失严重。所以使用金字塔前最好先把图像用padarray补成 2 的幂次倍或者只降两层不要贪多。4.2 顶层聚类结果向底层的映射策略从最小层开始聚类后要把中心传给下一层继续迭代。这个“用上层中心初始化下层”的过程比在下层随机初始化快得多而且不容易落入不同的局部最优。常见写法如下function [labels, centers] segmentWithPyramid(img, K, levels) pyr createPyramid(img, levels); nLevels length(pyr); centers []; % 初始无中心 labels []; for l nLevels:-1:1 feat preprocess(pyr{l}); if isempty(centers) % 第一次在最小层上完整迭代 [labels, centers] kmeansClusters(feat, K, 40); else % 后续层用上层中心初始化迭代次数可以降低 [labels, centers] kmeansClustersWithInit(feat, centers, 10); end end % 此时labels只对应pyr{1}的大小吗注意原图尺寸与pyr{1}一致 end注意pyr{1}是原图但preprocess后特征矩阵的行数等于原图宽乘高因此最终标签可以直接 reshape 回原图尺寸。这里容易犯错的是pyr{l}的尺寸每次缩小一半labels的尺寸也逐层变化需要保证在下一层初始化时把上一层中心的坐标映射到当前层而不是直接使用上一层的中心值。对于纯颜色特征中心值可以直接复用但如果特征里包含像素位置就需要按缩放比例把位置坐标乘 2 或除 2否则聚类中心会定位到错误的空间区域。4.3 K 值基准表与手肘法快速估计K 值是 k-means 彩色图像分割里最敏感的参数。下面是这套测试图里比较靠谱的参考区间测试图场景特点推荐 K 值注意事项starfish.jpg、beans.jpg单一前景背景颜色单一2-3K 太大背景会被拆碎redflowers2.jpg、flowers.jpg高饱和花朵与绿色叶3-5光照阴影区可能独立成簇buildings.jpg、angkorwat2.jpg天空、建筑、植被混合5-7天空与阴影常互相污染bagofgrain.jpg、patchofgrass.jpg纹理密集颜色相近7-9必须配合簇合并否则输出很碎如果不确定 K可以写一个最笨但有效的手肘法跑 K 从 2 到 10每次记录calculateQ.m输出的 Q 值或簇内误差平方和画曲线找拐点。拐点就是继续增大 K 时收益骤降的位置。注意每次跑之前固定随机种子否则曲线会毛刺很多无法判断拐点。4.4 金字塔层数和 K 的联动调参很多使用者把PyramidLevels当成“质量开关”以为层数越多越好。实际不是这样。金字塔层数增加的是空间平滑范围而不是分割精度。层数太多会让顶层图像丢失细长物体比如cornt.jpg里的玉米轮廓在降采样后可能只有两三个像素宽聚类中心根本不会为其单独建簇。一个可用的折中是PyramidLevels2加K比预期多 1然后用clusterMerge.m合并。这样既利用了金字塔提速又不会因为 K 太小漏掉细节区域。如果分割边界出现锯齿优先关掉金字塔而不是加层数直接在原图上跑 k-means再用imclose做一次形态学平滑。5. 簇合并、Q 值评估与三个立即可用的调参技巧5.1 clusterMerge.m 防止过分割的一次处理k-means 不感知空间位置纹理密集区域经常被拆成无数小簇。clusterMerge.m通过合并颜色中心距离极近的簇来减少这种现象。一个通用做法是计算簇中心两两距离不断合并最近且距离小于阈值的两个簇直到所有簇间距都大于阈值。合并时注意要用像素数加权平均而不是直接平均中心function labels mergeClusters(labels, centers, feat, threshold) % 合并距离过近的簇 K size(centers, 1); numPixels histcounts(labels, 1:K1); D squareform(pdist(centers)); D(1:K1:end) inf; while min(D(:)) threshold [c1, c2] find(D min(D(:)), 1); % 加权合并中心 n1 numPixels(c1); n2 numPixels(c2); newCenter (n1 * centers(c1, :) n2 * centers(c2, :)) / (n1 n2); % 更新标签 labels(labels c2) c1; centers(c1, :) newCenter; centers(c2, :) []; numPixels(c1) n1 n2; numPixels(c2) []; K K - 1; D squareform(pdist(centers)); D(1:K1:end) inf; end end阈值的选择没有万能公式。我一般先画合并前的分割区域图统计每个簇的像素占比。如果簇数量很多但最大簇只占 20%说明 K 设置过大合并阈值可以从最大簇中心距离的 10% 开始试。calculateQ.m中如果 Q 值在合并后不降反升说明这次合并是合理的。5.2 calculateQ.m 用来做参数网格搜索calculateQ.m计算的是簇内紧凑度和簇间分离度的综合指标。用它跑一组参数网格可以比人眼更客观地筛选参数组合。下面是一个可用脚本rng(42); Ks 2:7; threshes [0.05 0.1 0.2 0.4]; bestQ -inf; for K Ks for t threshes [labels, centers] segment_img(buildings.jpg, K, t); q computeQ(buildings.jpg, labels); fprintf(K%d thr%.2f Q%.4f\n, K, t, q); if q bestQ bestQ q; bestK K; bestT t; end end end实际使用中不要只看 Q 的最大值还要看排名靠前的几组参数对应的分割图。有些参数 Q 值很高但会把阴影全部并进天空区域视觉上几乎看不出建筑轮廓。这种情况下 Q 值就成了陷阱。5.3 固定随机种子与轮廓系数快速排查初始化问题最后一个立即可用的技巧是在所有入口脚本开头加上rng(42)。kmeansClusters.m用 randi 选初始中心不固定随机种子的话每次运行分割结果都不同。你在调 K 值时如果发现“脚本没改但两次输出不一样”基本都是随机种子的问题。排查初始化是否陷入局部最优用轮廓系数比用 Q 值更直接。matlab 自带silhouette函数输入特征矩阵和标签返回每个样本的轮廓值。轮廓值接近 1说明样本离邻近簇很遥远接近 0说明样本在两个簇的边界上。如果某张图多次运行时轮廓系数波动超过 0.2说明初始中心影响过大可以把maxIter提高到 60并多跑几次取 Q 值最高的那次结果。把上面的参数和代码在这套测试图上轮一遍你会看到从starfish.jpg的简单分割到patchofgrass.jpg的纹理分割差距主要在 K 值和合并阈值上而聚类算法本身只是提供了一个稳定的距离度量框架。这些调参经验迁移到自己的 matlab 彩色图像分割任务时可以先拿最小图验证再逐步放大避免在大图上反复试错浪费时间。本文还有配套的精品资源点击获取