MATLAB实现基于局部质心的无监督图像分割：2D/3D实战教程

在图像处理与分析领域，图像分割是一项基础且关键的任务，其目标是将图像划分为多个具有相似属性（如颜色、纹理、强度）的区域或对象。传统的分割方法往往依赖于大量的标注数据，这在医学影像、遥感等领域获取成本极高。因此，无监督图像分割方法，即不依赖人工标注即可自动发现图像结构的技术，具有重要的研究价值和实用意义。其中，基于局部质心的分割算法因其原理直观、实现相对简单，成为入门和探索无监督分割的一个经典切入点。

本文将以 MATLAB 为工具，深入讲解一种基于局部质心的无监督图像分割方法的核心原理，并提供完整的 2D 和 3D 图像分割实战代码。无论你是刚接触图像处理的同学，还是希望寻找一种轻量级分割方案的开发者，都能通过本文理解算法思想，并直接复现分割效果。

1. 图像分割与无监督学习核心概念

在深入代码之前，我们有必要厘清几个核心概念，这有助于理解我们即将实现的方法为何有效，以及它的适用边界。

1.1 什么是图像分割？

图像分割可以理解为“像素归类”问题。给定一张图像，分割算法需要为每一个像素分配一个唯一的标签，使得具有相同标签的像素在视觉上属于同一个物体或区域。例如，在一张医学CT图像中，分割的目标可能是将骨骼、软组织、背景等区分开来。

1.2 监督学习 vs. 无监督学习

这是机器学习的两大范式，在图像分割中体现得尤为明显：

• 监督学习分割：如 U-Net、Mask R-CNN 等深度学习方法。它们需要大量“图像-标注掩膜”配对的数据进行训练。模型学习从图像到分割结果的映射函数。优点是精度高，适用于定义明确、数据充足的场景（如“口腔疾病图像分割系统”）。
• 无监督学习分割：不需要任何标注数据。算法依据图像自身的统计特性、空间关系或像素相似性来发现“自然”的簇或边界。基于局部质心的方法、聚类算法（如K-means）、基于图的分割（如Graph Cut）等都属于此类。其优势在于无需标注，通用性强，但分割精度和语义一致性通常不如监督方法。

1.3 什么是“局部质心”？

“质心”通常指一个区域所有点的平均位置。在图像分割语境下，“局部质心”指的是：对于图像中的每一个像素，考察其周围一个局部邻域（例如一个 5x5 的窗口）内所有像素的特征（如灰度值、颜色向量），计算这个邻域内特征的平均值。这个平均值向量就被视为该像素所在“局部区域”的质心特征。

算法的核心思想：图像中属于同一物体的区域，其内部像素的局部邻域特征应该是相似的，因此它们的“局部质心”特征也会聚集在一起。相反，不同物体边界处的像素，其局部邻域会跨越不同区域，导致其“局部质心”特征与相邻像素差异较大。通过计算每个像素与其局部质心的差异，并以此作为分割依据，就能实现无监督的分割。

1.4 2D 与 3D 图像分割

• 2D 图像分割：处理的是二维矩阵，每个像素点由 (x, y) 坐标和强度值（或RGB向量）定义。这是最常见的形式。
• 3D 图像分割：处理的是三维体数据（如CT、MRI扫描结果），每个体素由 (x, y, z) 坐标和强度值定义。3D分割需要考虑体素之间的三维空间关系，计算复杂度更高，但对医学图像分析（“医学图像分割”、“3D点云”处理）至关重要。

我们实现的方法将同时支持 2D 和 3D 数据，其核心逻辑是相通的，只是邻域定义从二维窗口扩展到了三维立方体。

2. 环境准备与MATLAB基础

2.1 MATLAB 环境要求

本文代码基于 MATLAB R2018b 及以上版本编写和测试，主要使用了基本的矩阵运算和图像处理函数。确保你的 MATLAB 已安装以下工具箱（通常默认安装）：

• Image Processing Toolbox：用于图像读写、显示和基本处理。
• 无需其他特殊工具箱。

你可以通过以下命令检查：

% 检查Image Processing Toolbox是否安装 v = ver; if any(strcmp('Image Processing Toolbox', {v.Name})) disp('Image Processing Toolbox 已安装。'); else disp('警告：未找到Image Processing Toolbox，部分函数可能无法使用。'); end

2.2 项目文件结构建议

为了代码清晰，建议按如下结构组织你的项目文件夹：

local_centroid_segmentation/ ├── images/ % 存放测试图像 │ ├── test_2d.jpg │ └── test_3d.mat % 3D体数据通常保存为.mat文件 ├── utils/ % 工具函数 │ └── computeLocalCentroid.m ├── segment2D.m % 2D图像分割主函数 ├── segment3D.m % 3D图像分割主函数 ├── demo_2d.m % 2D演示脚本 └── demo_3d.m % 3D演示脚本

3. 算法原理与核心步骤拆解

基于局部质心的分割算法可以概括为以下四个核心步骤：

3.1 步骤一：图像预处理与特征提取

原始图像可能包含噪声，直接处理会影响质心计算的稳定性。通常先进行高斯滤波等平滑操作。对于彩色图像，可能需要将其转换为灰度图，或使用颜色特征。在本实现中，我们以灰度强度作为特征。

3.2 步骤二：计算局部邻域质心

这是算法的核心。对于图像中的每一个像素I(i, j)，我们定义一个大小为[ws, ws]（2D）或[ws, ws, ws]（3D）的滑动窗口。计算该窗口内所有像素特征值的平均值，作为该像素的局部质心C(i, j)。

% 伪代码逻辑 for i = 1:height for j = 1:width window = I(i-ws/2:i+ws/2, j-ws/2:j+ws/2); % 获取邻域 centroid(i, j) = mean(window(:)); % 计算均值（质心） end end

实际实现中，我们会使用imfilter或convn函数进行高效的卷积操作来替代循环，极大提升速度。

3.3 步骤三：构建差异图

计算原始图像I与局部质心图C的绝对差异D = |I - C|。在均匀区域，像素值与其局部质心接近，D值小；在边缘或纹理复杂区域，差异D值大。因此，差异图D可以看作是一个“边缘响应”或“非均匀性”的度量。

3.4 步骤四：基于差异图进行分割

得到差异图D后，有多种方式可以将其转换为分割结果：

1. 简单阈值法：设定一个阈值T，将D > T的像素标记为边界，D <= T的像素标记为内部。但这种方法只能得到二值边界，而非区域。
2. 聚类法（更常用）：将每个像素的特征表示为[I, D]或[I, C, D]，然后使用 K-means 等聚类算法对所有像素点进行聚类。属于同一簇的像素被赋予相同标签，实现分割。本文示例将采用这种方法。
3. 区域生长法：以差异较小的像素作为种子点，向周围相似区域生长。

4. 完整实战案例：2D 图像分割

我们将从一个具体的 2D 图像例子开始，逐步实现整个流程。

4.1 准备测试图像

你可以使用 MATLAB 自带的图像或任意你自己的图片。这里我们使用一张纹理图像。

% demo_2d.m clear; close all; clc; % 1. 读取或生成测试图像 % 使用内置图像 I = imread('cameraman.tif'); % 经典灰度图 % 或者使用纹理图像 % I = checkerboard(30, 4, 4); % 生成一个棋盘格图像 % I = uint8(255 * mat2gray(I)); figure(1); imshow(I); title('原始 2D 图像');

4.2 实现局部质心计算函数

我们将计算局部质心的功能封装成一个独立的函数，便于 2D 和 3D 复用。

% utils/computeLocalCentroid.m function centroid = computeLocalCentroid(image, windowSize) % 计算图像的局部质心图 % 输入： % image: 输入图像（2D或3D矩阵） % windowSize: 邻域窗口大小，标量（用于2D，如5）或向量（用于3D，如[5,5,5]） % 输出： % centroid: 与image同大小的局部质心图 % % 原理：使用均值滤波计算局部邻域的平均值。 % 创建均值滤波核 if isscalar(windowSize) % 2D 情况 h = fspecial('average', windowSize); centroid = imfilter(double(image), h, 'symmetric', 'same'); else % 3D 情况 % 创建一个三维的均值滤波核 kernel = ones(windowSize) / prod(windowSize); centroid = convn(double(image), kernel, 'same'); end end

4.3 2D 分割主函数

现在，我们编写主分割函数，集成所有步骤。

% segment2D.m function [labels, diffMap, centroidMap] = segment2D(I, windowSize, numClusters) % 基于局部质心的2D图像无监督分割 % 输入： % I: 输入灰度图像 (2D矩阵) % windowSize: 计算局部质心的邻域大小，必须为奇数，如 5, 7, 9 % numClusters: K-means聚类数目，即期望分割的区域数 % 输出： % labels: 分割标签图，大小与I相同 % diffMap: 原始图像与局部质心的差异图 % centroidMap: 局部质心图 % 步骤1: 转换为双精度浮点以便计算 I_double = double(I); % 步骤2: 计算局部质心图 centroidMap = computeLocalCentroid(I_double, windowSize); % 步骤3: 计算差异图 (绝对差异) diffMap = abs(I_double - centroidMap); % 步骤4: 特征构建与聚类分割 % 将每个像素的原始强度、局部质心和差异作为特征 [rows, cols] = size(I_double); features = [I_double(:), centroidMap(:), diffMap(:)]; % 每一行是一个像素的3维特征 % 使用K-means进行聚类 % ‘Replicates’参数设置多次随机初始聚类以避免局部最优，可根据需要调整 opts = statset('Display', 'final', 'MaxIter', 200); [labelIdx, ~] = kmeans(features, numClusters, 'Distance', 'sqeuclidean', ... 'Replicates', 3, 'Options', opts); % 将一维标签索引重塑为二维标签图 labels = reshape(labelIdx, rows, cols); % 可选：对标签图进行简单的形态学后处理，去除小噪声区域 % labels = medfilt2(labels, [3, 3]); end

4.4 运行与结果可视化

在演示脚本中调用主函数并展示结果。

% demo_2d.m (续) % 2. 设置算法参数 winSize = 7; % 局部邻域窗口大小，推荐奇数，如5,7,9。越大越平滑，但边界越模糊。 numClusters = 4; % 期望分割出的区域数量。需要根据图像内容先验估计。 % 3. 执行分割 [segLabels, diffMap, centroidMap] = segment2D(I, winSize, numClusters); % 4. 可视化结果 figure(2); subplot(2,2,1); imshow(I, []); title('原始图像'); subplot(2,2,2); imshow(centroidMap, []); title('局部质心图'); subplot(2,2,3); imshow(diffMap, []); title('差异图 |I-C|'); subplot(2,2,4); imagesc(segLabels); axis image; title('分割结果 (标签图)'); colormap(jet(numClusters)); colorbar; % 5. 将分割结果以彩色覆盖图形式显示在原图上 figure(3); imshow(I); hold on; % 生成一个随机的颜色映射给每个标签 randColors = rand(numClusters, 3); h = imagesc(segLabels); set(h, 'AlphaData', 0.4); % 设置部分透明度 colormap(randColors); title('分割区域叠加显示'); hold off; disp('2D 图像分割完成。');

4.5 结果分析与参数讨论

运行上述代码，你将看到四张图：

1. 原始图像：输入。
2. 局部质心图：比原图更平滑，细节被模糊，反映了每个像素邻域的平均亮度。
3. 差异图：高亮显示了原始图像与平滑后质心图的差异，通常对应边缘、纹理和噪声。
4. 分割标签图/叠加图：最终的分割结果，不同颜色代表算法识别出的不同区域。

关键参数影响：

• windowSize：控制局部邻域范围。值越小，质心图越能保留细节，差异图对噪声更敏感；值越大，平滑效果越强，可能模糊真实边界。通常尝试 5, 7, 9。
• numClusters：K-means 的簇数。这需要你对图像中感兴趣的区域数量有一个大致的先验估计。设置不当会导致“过分割”（区域太多）或“欠分割”（区域太少）。可以尝试使用“肘部法则”或轮廓系数来辅助选择，但对于无监督方法，这本身就是一个挑战。

5. 完整实战案例：3D 图像分割

3D 分割的逻辑与 2D 完全一致，只是数据维度增加，计算量更大。我们通常处理的是.mat文件或 DICOM 序列存储的体数据。

5.1 准备测试3D数据

由于公开3D图像数据不易获取，我们可以用 MATLAB 合成一个简单的3D体数据，包含几个不同强度的球体。

% demo_3d.m clear; close all; clc; % 1. 合成一个简单的3D体数据 (128x128x64) volSize = [128, 128, 64]; I_3d = zeros(volSize, 'uint8'); % 在体数据中创建几个不同强度的“球体” [x, y, z] = meshgrid(1:volSize(1), 1:volSize(2), 1:volSize(3)); center1 = [30, 30, 20]; radius1 = 15; center2 = [90, 90, 40]; radius2 = 20; center3 = [60, 60, 50]; radius3 = 10; sphere1 = sqrt((x-center1(1)).^2 + (y-center1(2)).^2 + (z-center1(3)).^2) <= radius1; sphere2 = sqrt((x-center2(1)).^2 + (y-center2(2)).^2 + (z-center2(3)).^2) <= radius2; sphere3 = sqrt((x-center3(1)).^2 + (y-center3(2)).^2 + (z-center3(3)).^2) <= radius3; I_3d(sphere1) = 150; % 中等灰度球体 I_3d(sphere2) = 50; % 暗色球体 I_3d(sphere3) = 220; % 亮色球体 % 添加一些高斯噪声模拟真实数据 I_3d = imnoise(I_3d, 'gaussian', 0, 0.01); disp(['3D体数据大小：', num2str(size(I_3d))]); % 显示中间切片 figure(1); imshow(I_3d(:,:,round(volSize(3)/2)), []); title('3D体数据中间切片 (XY平面)');

5.2 3D 分割主函数

3D函数与2D函数结构高度相似，主要区别在于卷积核是三维的，且特征矩阵的构建和重塑需要考虑三维。

% segment3D.m function [labels_3d, diffMap_3d, centroidMap_3d] = segment3D(V, windowSize3D, numClusters) % 基于局部质心的3D图像无监督分割 % 输入： % V: 输入3D体数据 (3D矩阵) % windowSize3D: 三维邻域大小，如 [5,5,5] 或标量5（表示[5,5,5]） % numClusters: K-means聚类数目 % 输出： % labels_3d: 3D分割标签体数据 % diffMap_3d: 3D差异体数据 % centroidMap_3d: 3D局部质心体数据 if isscalar(windowSize3D) windowSize3D = [windowSize3D, windowSize3D, windowSize3D]; end V_double = double(V); % 计算3D局部质心 centroidMap_3d = computeLocalCentroid(V_double, windowSize3D); % 计算3D差异图 diffMap_3d = abs(V_double - centroidMap_3d); % 构建特征矩阵 [dimX, dimY, dimZ] = size(V_double); numVoxels = dimX * dimY * dimZ; % 将三维特征展开成二维矩阵 (numVoxels x 3) features_3d = [V_double(:), centroidMap_3d(:), diffMap_3d(:)]; % 使用K-means聚类 % 注意：3D数据体素多，K-means计算可能很慢。可以考虑： % 1. 对体数据进行下采样后再聚类。 % 2. 使用更快的聚类算法，如MiniBatchKMeans (需要Statistics and Machine Learning Toolbox)。 % 3. 只使用部分体素样本进行聚类，然后插值。 disp('正在进行3D K-means聚类，数据量大时可能较慢...'); opts = statset('Display', 'final', 'MaxIter', 100); [labelIdx_3d, ~] = kmeans(features_3d, numClusters, 'Distance', 'sqeuclidean', ... 'Replicates', 2, 'Options', opts); % Replicates减少以加速 % 重塑标签 labels_3d = reshape(labelIdx_3d, dimX, dimY, dimZ); disp('3D分割完成。'); end

5.3 运行与3D结果可视化

3D结果的可视化比2D复杂，通常查看几个正交切片。

% demo_3d.m (续) % 2. 设置算法参数 winSize3D = 5; % 3D邻域，可以是标量5（表示[5,5,5]）或向量[5,5,5] numClusters3D = 4; % 我们合成了3个球体+背景，共4类 % 3. 执行3D分割 [segLabels3D, diffMap3D, centroidMap3D] = segment3D(I_3d, winSize3D, numClusters3D); % 4. 可视化结果 (显示中间切片) sliceZ = round(volSize(3)/2); sliceY = round(volSize(2)/2); sliceX = round(volSize(1)/2); figure(2); % 原始数据切片 subplot(2,3,1); imshow(I_3d(:,:,sliceZ), []); title('原始数据 (XY切片)'); subplot(2,3,2); imshow(squeeze(I_3d(:,sliceY,:)), []); title('原始数据 (XZ切片)'); subplot(2,3,3); imshow(squeeze(I_3d(sliceX,:,:)), []); title('原始数据 (YZ切片)'); % 分割结果切片 subplot(2,3,4); imagesc(segLabels3D(:,:,sliceZ)); axis image; title('分割标签 (XY)'); subplot(2,3,5); imagesc(squeeze(segLabels3D(:,sliceY,:))); axis image; title('分割标签 (XZ)'); subplot(2,3,6); imagesc(squeeze(segLabels3D(sliceX,:,:))); axis image; title('分割标签 (YZ)'); colormap(jet(numClusters3D)); % 5. 使用 isosurface 进行3D渲染展示（可选，更直观但计算稍慢） figure(3); % 为每个标签（除了背景，假设标签1是背景）绘制等值面 for k = 2:numClusters3D % 创建一个二值体，当前标签为1，其余为0 binaryVol = (segLabels3D == k); if any(binaryVol(:)) % 如果该标签存在 % 平滑一下等值面 binaryVol = smooth3(binaryVol, 'gaussian', 3); p = patch(isosurface(binaryVol, 0.5)); p.FaceColor = rand(1,3); p.EdgeColor = 'none'; p.FaceAlpha = 0.6; hold on; end end axis vis3d; grid on; view(3); camlight; lighting gouraud; title('3D分割结果等值面渲染'); xlabel('X'); ylabel('Y'); zlabel('Z'); hold off; disp('3D 图像分割完成。');

6. 常见问题与排查思路

在实际运行上述代码时，你可能会遇到一些典型问题。下表列出了常见问题及其解决方法：

7. 最佳实践与工程建议

基于局部质心的无监督分割方法简单有效，但在实际工程应用中，为了获得更鲁棒、更实用的结果，可以考虑以下优化方向和实践建议：

7.1 预处理至关重要

• 去噪：在计算局部质心前，应用适度的平滑滤波（如高斯滤波）可以显著抑制噪声对质心计算的影响，使差异图更能反映真实的结构边界，而非噪声点。
• 对比度增强：如果图像整体对比度较低，可以先进行直方图均衡化或对比度拉伸，增强区域间的差异，有助于后续聚类。
• 多尺度特征：单一尺度的windowSize可能无法同时捕捉大区域和小细节。可以尝试在多个尺度上计算局部质心和差异图，然后将这些多尺度特征融合，再送入聚类算法。

7.2 特征工程与聚类优化

• 特征选择：除了[I, C, D]，可以考虑加入纹理特征（如局部二值模式LBP、灰度共生矩阵特征）、梯度信息等，构建更丰富的特征向量，提升对复杂纹理的分割能力。
• 聚类算法选择：K-means 简单但需要指定K值且对噪声和初始值敏感。可以尝试：
  • 均值漂移 (Mean Shift)：无需指定聚类数量，能自动发现模态。
  • DBSCAN：基于密度，能发现任意形状的簇，并识别噪声点。
  • 谱聚类 (Spectral Clustering)：基于图论，在处理非凸数据分布时表现更好。
• 后处理：聚类得到的初始标签图可能存在小区域的孤立点或空洞。可以使用形态学操作（如开运算、闭运算）或连通组件分析来清理结果，合并过小的区域或填充孔洞。

7.3 针对3D数据的特殊优化

• 各向异性处理：医学影像等3D数据在Z轴（切片方向）的分辨率可能远低于XY平面。计算局部邻域时，应考虑各向异性，使用不同的窗口大小，例如[5,5,3]。
• 分块处理：对于超大的3D数据，可以将其分成重叠的小块分别处理，再合并结果，注意处理边界区域。
• 利用先验知识：在特定领域（如脑部MRI分割），可以结合解剖图谱等先验知识来约束聚类过程或解释聚类结果。

7.4 集成到完整流程

无监督分割结果通常可以作为更高级任务的起点：

• 监督学习的初始标注：为需要大量标注数据的深度学习模型提供初始的、粗糙的标注，再由人工进行精修，可以大幅减少人工标注工作量。
• 目标检测的候选区域：分割出的连通区域可以作为目标检测算法的候选框（Region Proposal）。
• 图像配准的预处理：分割出感兴趣区域后，可以只在该区域上进行图像配准，提高效率和精度。

本文提供的代码是一个完整的、可运行的起点。它清晰地展示了基于局部质心的无监督分割的核心流程。你可以以此为基线，根据具体的应用场景和图像特性，尝试上述的优化策略，逐步构建一个更强大、更鲁棒的分割工具。图像分割是一个广阔的领域，从传统的无监督方法到如今火热的深度学习，各有其适用场景。理解像本文这样的基础方法，对于深入掌握更复杂的模型有着不可替代的价值。