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部分函数可能无法使用。); end2.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:iws/2, j-ws/2:jws/2); % 获取邻域 centroid(i, j) mean(window(:)); % 计算均值质心 end end实际实现中我们会使用imfilter或convn函数进行高效的卷积操作来替代循环极大提升速度。3.3 步骤三构建差异图计算原始图像I与局部质心图C的绝对差异D |I - C|。在均匀区域像素值与其局部质心接近D值小在边缘或纹理复杂区域差异D值大。因此差异图D可以看作是一个“边缘响应”或“非均匀性”的度量。3.4 步骤四基于差异图进行分割得到差异图D后有多种方式可以将其转换为分割结果简单阈值法设定一个阈值T将D T的像素标记为边界D T的像素标记为内部。但这种方法只能得到二值边界而非区域。聚类法更常用将每个像素的特征表示为[I, D]或[I, C, D]然后使用 K-means 等聚类算法对所有像素点进行聚类。属于同一簇的像素被赋予相同标签实现分割。本文示例将采用这种方法。区域生长法以差异较小的像素作为种子点向周围相似区域生长。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 end4.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]); end4.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 结果分析与参数讨论运行上述代码你将看到四张图原始图像输入。局部质心图比原图更平滑细节被模糊反映了每个像素邻域的平均亮度。差异图高亮显示了原始图像与平滑后质心图的差异通常对应边缘、纹理和噪声。分割标签图/叠加图最终的分割结果不同颜色代表算法识别出的不同区域。关键参数影响windowSize控制局部邻域范围。值越小质心图越能保留细节差异图对噪声更敏感值越大平滑效果越强可能模糊真实边界。通常尝试 5, 7, 9。numClustersK-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分割完成。); end5.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. 常见问题与排查思路在实际运行上述代码时你可能会遇到一些典型问题。下表列出了常见问题及其解决方法问题现象可能原因解决思路MATLAB 报错未定义函数 ‘fspecial’Image Processing Toolbox 未安装。使用ver命令检查工具箱是否安装。如未安装需通过MATLAB附加功能管理器安装该工具箱。2D分割结果全是噪声没有连贯区域1.windowSize太小对噪声敏感。2.numClusters设置过大导致过分割。3. 图像本身噪声过大。1. 增大windowSize如从3改为7或9。2. 减小numClusters或尝试不同的值。3. 在计算质心前先对图像进行高斯滤波预处理 (imgaussfilt)。分割边界非常粗糙、呈块状windowSize太大导致局部质心过度平滑丢失了真实的细节边界。减小windowSize如从11改为5或7。需要在平滑噪声和保留边界之间权衡。3D分割代码运行极其缓慢甚至内存不足3D体数据体素数量巨大直接展开成特征矩阵进行K-means内存和计算开销都很大。1.降采样先对体数据V进行各向同性的降采样 (imresize3)在低分辨率上分割再将结果上采样回原尺寸。2.特征简化只使用[I, C]两维特征或仅使用I和D。3.采样聚类随机抽取一部分体素如10%进行K-means聚类得到聚类中心后为所有体素分配最近中心的标签。4. 使用MiniBatchKMeans(fitckmeans函数需Statistics and Machine Learning Toolbox)。K-means 结果每次运行都不一样K-means 算法对初始聚类中心敏感具有随机性。这是正常现象。通过增加‘Replicates’参数如从3增加到5或10让算法多次运行并选择最佳结果可以提高稳定性。但会增加计算时间。对于彩色图像如何应用上述代码默认处理灰度图像。将彩色图像转换为合适的颜色空间如Lab对每个通道分别计算局部质心或将RGB向量视为一个3维特征点计算其局部平均向量。这需要修改computeLocalCentroid函数以支持多通道输入。如何评价分割结果的好坏无监督分割缺乏真实标注Ground Truth定量评价困难。1.视觉评估观察分割区域是否与视觉感知一致。2.内部指标计算聚类内部的紧密度和类间的分离度如轮廓系数 (silhouette)。3.模拟数据在合成数据如本文的3D球体上可以计算与真实标签的吻合度如Dice系数。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。图像配准的预处理分割出感兴趣区域后可以只在该区域上进行图像配准提高效率和精度。本文提供的代码是一个完整的、可运行的起点。它清晰地展示了基于局部质心的无监督分割的核心流程。你可以以此为基线根据具体的应用场景和图像特性尝试上述的优化策略逐步构建一个更强大、更鲁棒的分割工具。图像分割是一个广阔的领域从传统的无监督方法到如今火热的深度学习各有其适用场景。理解像本文这样的基础方法对于深入掌握更复杂的模型有着不可替代的价值。