Matlab实战：基于K-means聚类算法实现图像分割

当你面对一张复杂的医学影像、卫星遥感图或工业检测照片，需要快速识别出其中的关键区域时，你会怎么做？手动勾勒边界不仅效率低下，而且主观性强。这时，“图像分割”技术就成了解决问题的关键。而在众多分割方法中，基于聚类算法的无监督分割因其无需预先标注数据而备受青睐。

今天要深入探讨的，正是利用Matlab平台，结合经典的K-means聚类算法，实现一套完整的图像分割方案。很多人以为K-means只是个简单的数据分组工具，但在图像处理领域，它却能以极低的计算成本，将图像像素按颜色、纹理或位置特征自动归类，从而勾勒出物体的轮廓。本文不仅提供可运行的源码，更会拆解其背后的设计逻辑、实现细节以及你一定会遇到的“坑”。

读完本文，你将能彻底理解如何将一幅彩色或灰度图像，通过K-means算法分割成具有清晰语义的区域，并掌握在Matlab中从零实现这一过程的全套技能。无论你是进行学术研究、完成课程设计，还是解决工程中的图像分析问题，这套方法都提供了一个坚实可靠的起点。

1. 图像分割：为什么K-means是一个“聪明”的起点？

在深入代码之前，我们必须先回答一个根本问题：图像分割领域有U-Net、Mask R-CNN等深度学习明星算法，为什么还要关注传统的K-means？

关键在于问题场景与资源约束。深度学习分割模型固然强大，但其代价是：

1. 数据依赖：需要大量精确标注的训练数据，标注成本极高。
2. 计算成本：训练和部署需要GPU等算力支持。
3. 模型复杂度：调参、训练流程复杂，不易理解和控制。

而K-means聚类分割的核心优势在于无监督和轻量级：

• 无需标注：算法直接分析图像像素本身的特征（如颜色），自动发现数据中的自然分组。
• 原理直观：核心思想（物以类聚）简单，易于实现和调试。
• 计算高效：对于中等尺寸图像，在CPU上即可快速得到结果。
• 结果可解释：分割出的区域直接对应特征空间中的簇，便于后续分析。

它特别适合以下场景：

• 预处理与粗分割：为更精细的算法（如深度学习）提供初始区域建议。
• 颜色主导的图像：例如分割蓝天、绿植、建筑；医学图像中不同组织（需在特定颜色空间下）。
• 快速原型验证：在资源有限或需要快速验证想法时。
• 教学与理解：是理解图像分割和聚类算法思想的绝佳案例。

因此，基于K-means的图像分割绝非过时技术，而是一个在特定条件下极其高效和实用的工具。它的价值在于用最小的代价解决一类明确的问题。

2. 核心原理：图像如何被“聚类”？

理解K-means用于图像分割，需要完成两次思维转换。

第一次转换：从图像到数据点我们看到的图像，在计算机中是一个数字矩阵。对于一幅M x N的彩色图像，每个像素点是一个包含RGB三个通道值的向量[R, G, B]。我们可以把整幅图像看作一个包含M*N个数据点的数据集，每个数据点位于一个三维（RGB）的特征空间中。灰度图像则是每个像素一个强度值，特征空间是一维的。

第二次转换：从聚类到分割K-means算法的任务是在这个特征空间中，找到K个中心点（簇心），并将所有像素点分配给距离最近的簇心。算法完成后，所有被分配到同一个簇心的像素点，就构成了图像中的一个“分割区域”。因为颜色相近的像素在特征空间中距离也近，自然会被聚到同一类。

算法流程简述：

1. 初始化：随机选择K个像素点的特征值作为初始簇心。
2. 分配阶段：遍历每个像素，计算其与所有簇心的距离（常用欧氏距离），将其分配给最近的簇心。
3. 更新阶段：重新计算每个簇所有像素特征的平均值，将该平均值作为新的簇心。
4. 迭代：重复步骤2和3，直到簇心的位置不再发生显著变化，或达到最大迭代次数。

最终，图像中每个像素都被赋予了一个“簇标签”（如1, 2, ..., K）。将相同标签的像素显示为同一种颜色，就得到了分割结果图。

3. 环境准备：你的Matlab需要什么？

在开始编写代码前，请确保你的Matlab环境已就绪。本项目对Matlab版本要求较为宽松，R2016a及之后的版本均可良好运行，因为它们都包含了必要的图像处理和基本数学工具箱。

核心工具箱检查：大部分功能基于Matlab核心函数，但为了图像读写和显示，请确认：

1. 打开Matlab，在命令行输入ver，查看已安装的产品列表。
2. 确保Image Processing Toolbox（图像处理工具箱）已安装。这是进行图像显示（imshow）、格式转换和某些高级操作所必需的。如果你的版本没有，部分基础功能可能受限，但核心聚类算法仍可运行。

项目文件结构建议：创建一个清晰的工作目录，例如Kmeans_Image_Segmentation，并在其中组织以下文件：

Kmeans_Image_Segmentation/ ├── images/ # 存放待分割的测试图像 │ ├── test1.jpg │ └── test2.png ├── kmeans_segmentation.m # 主函数文件 ├── visualize_results.m # 结果可视化脚本（可选） └── README.txt # 项目说明

良好的文件结构是避免路径错误的第一步。

4. 实战拆解：四步实现K-means图像分割

让我们将理论转化为可执行的步骤。整个过程可以分解为四个清晰的阶段。

4.1 第一步：图像读取与预处理

分割结果的质量很大程度上取决于输入数据的质量。原始图像可能包含噪声、光照不均或无关细节。

% 读取图像 original_img = imread('images/test1.jpg'); % 替换为你的图像路径 % 转换为双精度浮点型，便于后续计算 img_double = im2double(original_img); % (可选) 图像平滑去噪，使用高斯滤波 % 较小的滤波器尺寸(如3x3)和标准差(如0.5)可以平滑噪声同时保留边缘 filter_size = 3; sigma = 0.5; smoothed_img = imgaussfilt(img_double, sigma, 'FilterSize', filter_size); % 显示原始图像 figure(‘Position‘, [100, 100, 800, 400]); % 设置图形窗口位置和大小 subplot(1,2,1); imshow(original_img); title(‘原始图像‘);

关键点：im2double将像素值从0-255的整数转换到[0, 1]的浮点数区间，这对聚类算法的数值稳定性至关重要。滤波是可选项，对于噪声明显的图像推荐使用。

4.2 第二步：特征提取与数据重塑

这是连接图像和聚类算法的桥梁。我们需要从每个像素中提取特征，并将图像数据重塑为算法需要的n_samples x n_features矩阵。

% 方法1：仅使用颜色特征 (RGB) [M, N, C] = size(smoothed_img); % 将图像重塑为 (M*N) 行， C 列的矩阵，每一行是一个像素的特征向量 color_features = reshape(smoothed_img, M*N, C); % 方法2：结合颜色和空间位置特征 (更常用，能产生空间连续的区域) [X, Y] = meshgrid(1:N, 1:M); % 生成像素坐标网格 % 将坐标归一化到 [0,1] 区间，使其与颜色特征量级匹配 X_norm = X / max(X(:)); Y_norm = Y / max(Y(:)); % 提取颜色特征 color_features = reshape(smoothed_img, M*N, C); % 将坐标特征与颜色特征拼接 spatial_features = [X_norm(:), Y_norm(:)]; % 加权拼接：颜色权重通常更高，以颜色为主导，空间信息辅助平滑 weight_color = 0.8; weight_spatial = 0.2; all_features = [weight_color * color_features, weight_spatial * spatial_features]; % 本次演示使用结合特征的方案 features_for_clustering = all_features;

深度解析：纯颜色特征可能导致空间上不连续但颜色相似的区域被分为一类。加入归一化的坐标特征 (X_norm,Y_norm) 后，算法会倾向于将空间位置相近的像素聚在一起，从而得到更紧凑、连续的分割块。权重参数 (weight_color,weight_spatial) 是重要的调节旋钮，需要根据图像内容调整。

4.3 第三步：应用K-means聚类

这是核心步骤。我们将使用Matlab内置的kmeans函数。选择K值（分割区域数）是本步骤最大的挑战。

% 设定聚类数目K K = 4; % 例如，想把图像分割成4个主要区域 % 设置K-means参数 max_iterations = 100; % 最大迭代次数 replicates = 3; % 重复聚类次数，使用最佳结果，避免局部最优 options = statset(‘MaxIter‘, max_iterations); % 执行K-means聚类 % ‘Display‘, ‘final‘ 显示最终收敛信息 % ‘Replicates‘, replicates 重复运行取最好结果 [cluster_idx, cluster_center] = kmeans(features_for_clustering, K, ... ‘Options‘, options, ... ‘Display‘, ‘final‘, ... ‘Replicates‘, replicates, ... ‘Distance‘, ‘sqeuclidean‘); % 使用平方欧氏距离 % cluster_idx: 一个 (M*N) x 1 的向量，包含每个像素点的簇标签(1到K) % cluster_center: 一个 K x n_features 的矩阵，每一行是一个簇的中心

参数精讲：

• K值选择：没有绝对标准。可以尝试使用“肘部法则”（Elbow Method），绘制不同K值对应的聚类误差（kmeans函数返回的到簇心距离之和），选择误差下降速度骤缓的点。对于简单图像，K=3或4通常是个好的开始。
• Replicates：由于K-means初始簇心随机，可能收敛到局部最优。设置replicates > 1让算法多次运行并选择误差最小的结果，增加稳定性。
• Distance：对于颜色和坐标这类连续数值特征，‘sqeuclidean‘（平方欧氏距离）是最常用且有效的选择。

4.4 第四步：标签重构与结果显示

将一维的簇标签向量重新变回二维图像格式，并可视化分割结果。

% 将簇标签向量重塑回与原图像相同大小的矩阵 segmented_labels = reshape(cluster_idx, M, N); % 创建一个彩色分割显示图像 % 为每个簇标签分配一个鲜明的颜色 colors = lines(K); % lines是Matlab的一个颜色图，生成K种区分度好的颜色 segmented_rgb = zeros(M, N, 3); for k = 1:K color_mask = (segmented_labels == k); for ch = 1:3 channel = segmented_rgb(:, :, ch); channel(color_mask) = colors(k, ch); segmented_rgb(:, :, ch) = channel; end end % 显示分割结果 subplot(1,2,2); imshow(segmented_rgb); title([‘K-means分割结果 (K=‘, num2str(K), ‘)‘]); % 可选：将原图与分割结果叠加显示，以观察边界 figure; imshow(original_img); hold on; % 只显示区域边界 boundaries = bwperim(segmented_labels ~= 0); % 获取标签区域的边界 visboundaries(boundaries, ‘Color‘, ‘r‘, ‘LineWidth‘, 1.5); title(‘原图叠加分割边界‘); hold off;

结果解读：第一幅结果图用不同色块直接显示了分割区域。第二幅叠加图在原图上用红线勾勒了分割边界，能更直观地评估分割精度，看边界是否与物体边缘吻合。

5. 完整源码与模块化函数设计

将上述步骤封装成一个健壮的、可复用的函数，是工程化的关键。下面提供一个更完整、更模块化的主函数示例。

% 文件：kmeans_image_segmentation.m function [segmented_labels, cluster_centers] = kmeans_image_segmentation(img_path, K, use_spatial, weight_color, spatial_weight) % 基于K-means的图像分割主函数 % 输入： % img_path - 字符串，输入图像路径 % K - 整数，聚类数目 % use_spatial - 逻辑值，true表示使用空间特征，false则仅用颜色 % weight_color - 标量，颜色特征的权重 (默认0.8) % spatial_weight- 标量，空间特征的权重 (默认0.2)，仅当use_spatial=true时有效 % 输出： % segmented_labels - MxN矩阵，每个像素的簇标签(1-K) % cluster_centers - K x n_features矩阵，聚类中心 % 1. 读取与预处理图像 if ~exist(img_path, ‘file‘) error(‘图像文件不存在: %s‘, img_path); end original_img = imread(img_path); img = im2double(original_img); [M, N, C] = size(img); % 可选：应用高斯滤波去噪（根据图像噪声情况开启） % img = imgaussfilt(img, 0.5); % 2. 特征提取 color_feat = reshape(img, M*N, C); if use_spatial % 生成并归一化空间坐标 [X, Y] = meshgrid(1:N, 1:M); X_norm = X / max(X(:)); Y_norm = Y / max(Y(:)); spatial_feat = [X_norm(:), Y_norm(:)]; % 特征加权拼接 % 确保权重和为1，或根据经验调整 if nargin < 4 weight_color = 0.8; spatial_weight = 0.2; end all_features = [weight_color * color_feat, spatial_weight * spatial_feat]; else all_features = color_feat; end % 3. K-means聚类 max_iter = 100; replicates = 5; % 增加重复次数以获得更稳定结果 opts = statset(‘MaxIter‘, max_iter); fprintf(‘开始K-means聚类，K=%d, 使用空间特征: %d...\n‘, K, use_spatial); [idx, centers] = kmeans(all_features, K, ‘Options‘, opts, ... ‘Replicates‘, replicates, ... ‘Display‘, ‘off‘); % 关闭迭代显示，保持整洁 % 4. 重塑标签矩阵 segmented_labels = reshape(idx, M, N); cluster_centers = centers; fprintf(‘聚类完成。\n‘); end

调用示例：

% 主脚本：main_demo.m clear; clc; close all; % 设置参数 image_file = ‘images/your_image.jpg‘; % 请替换为实际图像 K = 4; use_spatial_feature = true; color_w = 0.7; spatial_w = 0.3; % 执行分割 [labels, centers] = kmeans_image_segmentation(image_file, K, use_spatial_feature, color_w, spatial_w); % 可视化结果 original_img = imread(image_file); visualize_segmentation(original_img, labels, K);

可视化辅助函数：

% 文件：visualize_segmentation.m function visualize_segmentation(original_img, label_map, K) figure(‘Position‘, [50, 50, 1200, 500]); % 子图1：原始图像 subplot(1,3,1); imshow(original_img); title(‘原始图像‘, ‘FontSize‘, 12); % 子图2：伪彩色分割图 subplot(1,3,2); % 使用jet或hsv颜色图，区分度更好 cmap = jet(K); % 或 hsv(K), lines(K) rgb_label = label2rgb(label_map, cmap, ‘k‘, ‘shuffle‘); % ‘k‘背景为黑色 imshow(rgb_label); title([‘伪彩色分割图 (K=‘, num2str(K), ‘)‘], ‘FontSize‘, 12); % 子图3：原图叠加边界 subplot(1,3,3); imshow(original_img); hold on; % 为每个区域绘制边界 for k = 1:K region_boundary = bwperim(label_map == k); [B_y, B_x] = find(region_boundary); plot(B_x, B_y, ‘.‘, ‘Color‘, cmap(k,:), ‘MarkerSize‘, 1); % 用小点绘制边界 end hold off; title(‘原图叠加分割边界‘, ‘FontSize‘, 12); end

这个模块化设计将核心算法、参数配置和结果展示分离，使得代码更清晰、易于调试和扩展。

6. 运行效果与评估：如何判断分割得好不好？

运行上述代码后，你将得到类似下图的结果。但如何客观评价分割质量呢？

定性评估（肉眼观察）：

1. 区域一致性：同一物体内部是否被分割为同一个区域？有无过度分割（一个物体被分成多块）或欠分割（多个物体被合并成一块）？
2. 边界贴合度：分割边界是否与物体的真实边缘基本吻合？
3. 语义相关性：分割出的区域是否具有明确的语义（如“天空”、“道路”、“车辆”）？

定量评估（如有Ground Truth）：如果你有标准的分割标注图（真值），可以使用以下指标，这些也是学术论文中常用的：

% 假设 gt_labels 是标准真值标签图，pred_labels 是我们的预测标签图 % 注意：K-means产生的标签编号是任意的，需要与真值标签对齐后才能计算。 % 这里演示调整兰德指数(ARI)的计算（需要Statistics and Machine Learning Toolbox） % 首先将标签图展平为向量 gt_vector = double(gt_labels(:)); pred_vector = double(pred_labels(:)); % 计算调整兰德指数 (Adjusted Rand Index, ARI) % ARI取值范围[-1,1]，越大越好，1表示完全一致，0表示随机划分。 ari = rand_index(gt_vector, pred_vector, ‘adjusted‘); % 需要自定义或使用第三方函数 % 对于像素级精度 accuracy = sum(gt_labels(:) == aligned_pred_labels(:)) / numel(gt_labels); fprintf(‘像素精度: %.4f\n‘, accuracy);

对于无监督的K-means，更常用的内部评估指标是簇内距离和（WCSS），即kmeans函数可返回的每个点到其簇心距离的总和。绘制WCSS随K值变化的“肘部曲线”可以帮助选择K。

7. 常见问题与调试指南

在实际操作中，你几乎一定会遇到以下问题。这里提供系统的排查思路。

高级调试技巧：特征可视化在特征提取后、聚类前，可以抽样可视化特征空间，直观理解算法“看到”的数据。

% 抽样一部分像素点进行可视化 sample_indices = randperm(size(all_features, 1), 2000); % 随机取2000个点 sample_features = all_features(sample_indices, :); % 如果是RGB特征（3维），可以用3D散点图 if size(sample_features,2) == 3 && ~use_spatial figure; scatter3(sample_features(:,1), sample_features(:,2), sample_features(:,3), 10, ... img_double(sample_indices), ‘filled‘); % 用原色着色 xlabel(‘Red‘); ylabel(‘Green‘); zlabel(‘Blue‘); title(‘像素点在RGB特征空间中的分布‘); colorbar; end % 更高维特征可使用PCA降维后可视化

8. 超越基础：最佳实践与进阶思路

掌握了基础实现后，以下实践和思路能让你的分割项目更上一层楼。

8.1 特征工程：从RGB到更优的颜色空间

RGB颜色空间对光照变化敏感。HSV和CIELab是更常用的颜色空间。

% 转换为HSV颜色空间 img_hsv = rgb2hsv(original_img); % 通常使用H(色调)和S(饱和度)通道，丢弃V(明度)通道或谨慎使用 H = img_hsv(:,:,1); S = img_hsv(:,:,2); % 将H和S通道作为特征 hsv_features = [H(:), S(:)]; % 注意：H通道是角度（0-1代表0-360度），计算距离时需处理周期性问题 % 一种简单方法是计算角度差的最小值：min(|h1-h2|, 1-|h1-h2|)

8.2 自动确定K值：肘部法则的实现

通过遍历一系列K值，计算WCSS，找到拐点。

max_K_to_try = 10; wcss = zeros(max_K_to_try, 1); % 簇内距离平方和 for k = 1:max_K_to_try [~, ~, sumd] = kmeans(features_for_clustering, k, ‘Options‘, opts, ... ‘Replicates‘, 3, ‘Display‘, ‘off‘); wcss(k) = sum(sumd); % 所有簇的距离和 end % 绘制肘部曲线 figure; plot(1:max_K_to_try, wcss, ‘bo-‘); xlabel(‘聚类数目 K‘); ylabel(‘簇内距离平方和 (WCSS)‘); title(‘肘部法则曲线‘); grid on; % 观察曲线，选择WCSS下降速度突然变缓的点作为K

8.3 后处理：优化分割结果

K-means产生的原始标签图可能存在小噪声区域。可以使用形态学操作进行后处理。

% 假设 segmented_labels 是初步结果 for k = 1:K % 获取当前簇的二进制掩膜 mask = (segmented_labels == k); % 1. 开运算：先腐蚀后膨胀，去除小噪点 se = strel(‘disk‘, 2); % 创建一个半径为2的圆盘结构元素 mask_cleaned = imopen(mask, se); % 2. 闭运算：先膨胀后腐蚀，填充小孔洞 mask_cleaned = imclose(mask_cleaned, se); % 3. 移除小面积区域（可选） mask_cleaned = bwareaopen(mask_cleaned, 50); % 移除面积小于50像素的区域 % 更新标签 segmented_labels(mask_cleaned) = k; % 注意：被移除的小区域需要被重新分配标签，这里简化处理 end % 注意：后处理可能改变区域像素数，需谨慎评估。

8.4 工程化建议

• 批量处理：将主函数封装好后，使用循环或parfor并行循环处理一个文件夹下的所有图像。
• 参数调优：对于特定类型图像（如所有卫星图），可以设计一个小的验证集，用网格搜索法寻找最优的K、weight_color、spatial_weight等参数组合。
• 结果保存：不仅保存分割后的彩色图，也保存标签矩阵 (segmented_labels)，便于后续定量分析或作为其他算法的输入。
• 性能记录：在日志中记录每张图像的处理时间、使用的参数和最终WCSS，用于分析和比较。

从加载图像到得到分割结果，你已走完了基于K-means的图像分割全流程。这套方法的精髓在于将图像像素映射为特征空间中的点，再通过聚类发现其内在结构。它快速、直观且无需训练数据，是解决许多实际分割问题的利器。

然而，务必清楚其边界：K-means基于全局特征分布，对复杂纹理、遮挡物体和强度不均匀的目标分割能力有限。当面对这类挑战时，今天的实践成为了理解更高级算法（如基于图割的GrabCut、基于超像素的SLIC，乃至深度学习分割模型）的基石。你可以尝试将本文的K-means结果作为这些算法的初始化输入，探索传统与前沿方法的结合。

建议你将本文的源码作为工具箱保存，下次遇到需要从图像中快速提取颜色一致区域的任务时，它很可能就是最高效的解决方案。