如果你正在处理医学影像、遥感图像或任何缺乏标注数据的图像分割任务，那么“无监督”这三个字可能就是你的救星。传统的深度学习分割方法，如U-Net，虽然效果出色，但严重依赖大量、高质量的标注数据。获取这些数据不仅成本高昂，在医学、工业检测等领域甚至可能涉及隐私或专业壁垒。那么，有没有一种方法，能让计算机像人眼一样，仅凭图像自身的纹理、亮度和结构信息，就自动找出其中的物体边界？

答案是肯定的，而且其核心思想可能比你想象的更“物理”。本文要深入探讨的，正是这样一种基于**局部质心（Local Centroid）**的无监督图像分割方法。它不依赖任何预训练模型或标注数据，仅通过计算图像局部区域的“质量中心”来驱动像素的聚合与分离，巧妙地实现了对2D和3D图像的自动分割。

很多人初次接触“无监督分割”时，可能会联想到复杂的聚类算法（如K-means）或传统的边缘检测。但基于局部质心的方法提供了一个截然不同的视角：它将图像视为一个物理场，像素的灰度值视为“质量”，通过迭代计算每个像素所属的局部质心，让相似的像素“自然”地聚集到同一个质心周围，最终形成分割区域。这种方法在处理噪声图像、灰度不均匀图像时，往往表现出比简单阈值法或聚类法更强的鲁棒性。

本文将为你彻底拆解这个算法的Matlab实现。你不仅能获得可直接运行的完整代码，更能理解其背后的数学原理和迭代动力。我们会从最基础的概念讲起，一步步搭建环境、解析核心函数、运行示例，并深入探讨参数调优、常见问题与工程实践。无论你是医学图像分析的研究者，还是从事计算机视觉的工程师，这篇文章都将为你提供一个强大且无需标注数据的工具箱。

1. 核心问题：为什么需要无监督的局部质心分割？

在深入代码之前，我们必须先厘清一个根本问题：既然有那么多现成的分割工具，为什么还要关注这个相对“古老”的算法？

关键在于适用场景的不可替代性。考虑以下情况：

1. 标注数据稀缺或无法获取：许多特殊的医学影像（如罕见病）、工业缺陷样本、天文观测图像，根本没有现成的标注数据集。
2. 数据分布多变：如果你要处理的数据来源不一（如不同医院、不同设备拍摄的CT），其对比度、亮度差异很大，一个基于固定数据训练的监督模型可能泛化能力很差。
3. 需要快速原型验证：在项目初期，你需要一个快速、无需训练的分割工具来验证想法的可行性，而不是花几周时间去标注数据、训练模型。
4. 对可解释性要求高：局部质心法原理基于直观的物理类比（质心吸引），每个步骤的结果都可以可视化，便于调试和理解，不像深度神经网络那样是个“黑箱”。

基于局部质心的分割方法，正是为解决这些痛点而生。它放弃了从数据中“学习”模式，转而利用图像固有的空间与灰度一致性，通过迭代优化找到一个稳定的分割状态。它的优势不在于超越SOTA的精度，而在于零样本启动、强鲁棒性和高可解释性。

2. 算法原理：将图像视为一个物理系统

理解这个算法，我们可以做一个生动的类比：把一张灰度图像想象成一个凹凸不平的金属板，图像中每个像素点的灰度值（亮度）代表该点金属板的“高度”或“密度”。

• 像素= 金属板上的一个质点。
• 灰度值= 该质点的“质量”（越亮，质量越大）。
• 局部区域= 以某个质点为中心的一个小窗口（比如3x3, 5x5）。
• 局部质心= 在这个小窗口内，所有质点的“质量中心”。质量大的点（亮区）会把质心拉向自己。

算法的核心思想是：每个像素点都应该属于它所在局部区域的质心所代表的“物体”。通过迭代，每个像素点不断向它局部窗口计算出的质心位置“靠拢”，同时质心的位置又根据归属于它的像素点重新计算。最终，相似的像素点会收敛到同一个或几个稳定的质心，从而完成分割。

数学描述（以2D图像为例）： 对于一个像素点(i, j)，定义一个以其为中心的局部窗口W。该窗口内所有像素的坐标和灰度值的加权平均，即为局部质心(C_i, C_j)：

[ C_i = \frac{\sum_{(m,n) \in W} I(m, n) \cdot m}{\sum_{(m,n) \in W} I(m, n)}, \quad C_j = \frac{\sum_{(m,n) \in W} I(m, n) \cdot n}{\sum_{(m,n) \in W} I(m, n)} ]

其中，I(m, n)是像素点(m, n)的灰度值。

迭代过程：

1. 初始化：可以将每个像素的坐标作为其初始“标签”或“所属质心位置”。
2. 质心计算：对于当前每个像素，根据其当前所属的质心标签，找到所有共享同一标签的像素，在这些像素构成的区域上重新计算质心位置。
3. 标签更新：对于每个像素，在其周围局部窗口内，寻找距离哪个计算出的质心最近（或根据灰度相似性判断），就将自己的标签更新为该质心的标签。
4. 收敛判断：重复步骤2和3，直到所有像素的标签不再变化，或变化非常小，迭代停止。

这个过程类似于均值漂移（Mean Shift）聚类，但它的操作定义在规则的图像网格上，并且与区域的物理质心概念紧密结合。

扩展到3D：对于3D图像（如CT、MRI序列），原理完全一致，只是像素点变为体素（Voxel），局部窗口从2D正方形变为3D立方体，坐标(i, j)变为(i, j, k)。计算3D质心的公式是类似的加权平均。

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

本算法的实现和实验完全依赖于Matlab。以下是确保你能顺利运行代码的清单：

1. Matlab版本：建议使用R2018a及以上版本。算法主要使用基础矩阵运算和图像处理工具箱，这些功能在较旧的版本中也存在，但新版在内存管理和图形显示上更优。本文代码已在 R2021b 上测试通过。
2. 必需工具箱：
   • Image Processing Toolbox：这是核心，用于图像的读写（imread,imwrite）、显示（imshow）和基本的图像操作。
   • 不需要深度学习、计算机视觉或统计工具箱。这体现了该算法的轻量性。
3. 硬件：无特殊要求。对于大型3D图像（如256x256x100），迭代算法可能会消耗较多内存和计算时间。确保你有足够的RAM（建议8GB以上）。
4. 代码获取：你需要将本文后续章节提供的核心函数代码保存为.m文件。

项目文件结构建议： 创建一个项目文件夹，例如LocalCentroidSegmentation，内部结构如下：

LocalCentroidSegmentation/ ├── main_2d_demo.m % 2D图像分割主演示脚本 ├── main_3d_demo.m % 3D图像分割主演示脚本 ├── local_centroid_seg_2d.m % 2D分割核心函数 ├── local_centroid_seg_3d.m % 3D分割核心函数 ├── data/ │ ├── test_2d_image.png % 示例2D图像 │ └── test_3d_volume.mat % 示例3D体数据（.mat格式） └── results/ % 保存输出结果的文件夹

4. 核心函数实现与逐行解析

我们将分别实现2D和3D版本的核心函数。理解这些代码是掌握该方法的关键。

4.1 2D图像分割核心函数

将以下代码保存为local_centroid_seg_2d.m。

function [labels, centroids] = local_centroid_seg_2d(I, window_size, max_iters, convergence_thresh) % 基于局部质心的无监督2D图像分割 % 输入： % I - 输入灰度图像 (2D矩阵, double类型，范围建议[0,1]) % window_size - 局部窗口的半径（例如，3表示7x7的窗口，因为 2*3+1=7） % max_iters - 最大迭代次数 % convergence_thresh - 标签变化的收敛阈值（比例） % 输出： % labels - 分割标签图，大小与I相同，整数类型 % centroids - 最终迭代后每个标签对应的质心坐标列表 [x, y] [height, width] = size(I); num_pixels = height * width; % 1. 初始化：每个像素的标签为其自身的线性索引 labels = reshape(1:num_pixels, height, width); % 2. 将图像灰度值归一化并作为权重（质量） weights = (I - min(I(:))) / (max(I(:)) - min(I(:)) + eps); % 归一化到[0,1] weights = weights + 0.01; % 避免零权重，保证数值稳定性 for iter = 1:max_iters fprintf('迭代 %d/%d...\n', iter, max_iters); labels_old = labels; % 3. 计算每个唯一标签对应的质心 unique_labels = unique(labels(:)); centroids = zeros(length(unique_labels), 2); % 存储[x, y]坐标 for idx = 1:length(unique_labels) label = unique_labels(idx); mask = (labels == label); % 属于当前标签的像素掩码 % 计算该区域内的加权质心 [y_coords, x_coords] = find(mask); % 注意：find返回的是(row, col)，即(y, x) pixel_weights = weights(mask); if sum(pixel_weights) > 0 cent_x = sum(x_coords .* pixel_weights) / sum(pixel_weights); cent_y = sum(y_coords .* pixel_weights) / sum(pixel_weights); centroids(idx, :) = [cent_x, cent_y]; else centroids(idx, :) = [mean(x_coords), mean(y_coords)]; end end % 4. 为每个像素更新标签：在局部窗口内寻找最近的质心 new_labels = zeros(size(labels)); % 为了效率，我们遍历每个质心，计算其影响范围 % 更精确但低效的做法是遍历每个像素，在其窗口内搜索质心。这里采用近似高效方法。 % 创建坐标网格 [X, Y] = meshgrid(1:width, 1:height); for idx = 1:size(centroids, 1) cent_x = centroids(idx, 1); cent_y = centroids(idx, 2); label = unique_labels(idx); % 计算每个像素到该质心的距离 dist_map = sqrt((X - cent_x).^2 + (Y - cent_y).^2); % 创建一个窗口掩码：只考虑质心周围 window_size 范围内的像素 window_mask = (abs(Y - cent_y) <= window_size) & (abs(X - cent_x) <= window_size); % 在窗口内，如果当前像素距离这个质心比之前记录的更近，则更新标签 % 初始化时，new_labels为0，距离为无穷大 if iter == 1 current_min_dist = inf(size(labels)); else % 对于已在窗口内的像素，与当前最小距离比较 update_mask = window_mask & (dist_map < current_min_dist); new_labels(update_mask) = label; current_min_dist(update_mask) = dist_map(update_mask); end end % 处理未被任何窗口覆盖的像素（理论上很少，保留原标签） uncovered_pixels = (new_labels == 0); new_labels(uncovered_pixels) = labels_old(uncovered_pixels); labels = new_labels; % 5. 检查收敛：计算标签变化的像素比例 change_ratio = sum(labels(:) ~= labels_old(:)) / num_pixels; fprintf(' 标签变化比例: %.4f\n', change_ratio); if change_ratio < convergence_thresh fprintf('算法在 %d 次迭代后收敛。\n', iter); break; end end % 6. 重新映射标签为连续的整数（1,2,3,...），便于显示和后续处理 [labels, ~] = cleanup_labels(labels); end function [clean_labels, map] = cleanup_labels(labels) % 辅助函数：将标签重新映射为连续的整数 unique_vals = unique(labels(:)); map = containers.Map(unique_vals, 1:length(unique_vals)); clean_labels = zeros(size(labels)); for i = 1:numel(labels) clean_labels(i) = map(labels(i)); end end

关键逻辑解析：

• 初始化：每个像素独立为一个类，这是最细粒度的起点。
• 质心计算：在每次迭代中，算法根据当前的分割结果（labels），计算每个连通区域的灰度加权质心。weights由图像归一化得到，亮区像素拥有更高的“质量”，对质心位置有更大的拉力。
• 标签更新：这是算法的核心。我们采用了一种“由质心辐射”的策略。对于每个质心，我们只考虑其周围window_size范围内的像素。在这个局部窗口内，像素点会选择距离最近的质心作为自己的新标签。这比遍历每个像素再搜索全局质心要高效得多，也符合“局部”的原则。
• 收敛判断：当重新分配标签的像素比例低于阈值时，认为系统达到稳定状态。
• 标签清理：最终输出的标签可能是离散的整数（如1, 5, 100），cleanup_labels函数将其重映射为连续的（1, 2, 3, ...），方便可视化。

4.2 3D图像分割核心函数

将以下代码保存为local_centroid_seg_3d.m。其逻辑与2D版本完全一致，只是扩展到了三维。

function [labels, centroids] = local_centroid_seg_3d(V, window_size, max_iters, convergence_thresh) % 基于局部质心的无监督3D图像分割 % 输入： % V - 输入3D灰度体数据 (3D矩阵, double类型，范围建议[0,1]) % window_size - 局部窗口的半径（在3D中，这是立方体窗口的半径） % max_iters - 最大迭代次数 % convergence_thresh - 标签变化的收敛阈值（比例） % 输出： % labels - 分割标签体，大小与V相同，整数类型 % centroids - 最终迭代后每个标签对应的质心坐标列表 [x, y, z] [depth, height, width] = size(V); % 注意：Matlab中3D矩阵是 (行, 列, 页) -> (y, x, z)? 这里按常规(x,y,z)理解。 % 更常见的顺序：size(V) 返回 [height, width, depth]。我们需要确认。 % 假设输入 V 是 [Height, Width, Depth] [height, width, depth] = size(V); num_voxels = height * width * depth; labels = reshape(1:num_voxels, height, width, depth); weights = (V - min(V(:))) / (max(V(:)) - min(V(:)) + eps); weights = weights + 0.01; % 创建3D坐标网格 [X, Y, Z] = meshgrid(1:width, 1:height, 1:depth); for iter = 1:max_iters fprintf('3D迭代 %d/%d...\n', iter, max_iters); labels_old = labels; unique_labels = unique(labels(:)); centroids = zeros(length(unique_labels), 3); % 存储[x, y, z] % 计算每个标签的质心 for idx = 1:length(unique_labels) label = unique_labels(idx); mask = (labels == label); % 找到属于该标签的所有体素坐标 [y_coords, x_coords, z_coords] = ind2sub([height, width, depth], find(mask)); voxel_weights = weights(mask); if sum(voxel_weights) > 0 cent_x = sum(x_coords .* voxel_weights) / sum(voxel_weights); cent_y = sum(y_coords .* voxel_weights) / sum(voxel_weights); cent_z = sum(z_coords .* voxel_weights) / sum(voxel_weights); centroids(idx, :) = [cent_x, cent_y, cent_z]; else centroids(idx, :) = [mean(x_coords), mean(y_coords), mean(z_coords)]; end end % 更新标签：每个体素选择其局部窗口内最近的质心 new_labels = zeros(size(labels)); current_min_dist = inf(size(labels)); for idx = 1:size(centroids, 1) cent_x = centroids(idx, 1); cent_y = centroids(idx, 2); cent_z = centroids(idx, 3); label = unique_labels(idx); % 计算3D欧氏距离 dist_map = sqrt((X - cent_x).^2 + (Y - cent_y).^2 + (Z - cent_z).^2); % 创建3D局部窗口掩码 window_mask = (abs(Y - cent_y) <= window_size) & ... (abs(X - cent_x) <= window_size) & ... (abs(Z - cent_z) <= window_size); % 在窗口内，找到距离更近的体素并更新其标签 update_mask = window_mask & (dist_map < current_min_dist); new_labels(update_mask) = label; current_min_dist(update_mask) = dist_map(update_mask); end % 处理未被覆盖的体素 uncovered_voxels = (new_labels == 0); new_labels(uncovered_voxels) = labels_old(uncovered_voxels); labels = new_labels; % 收敛判断 change_ratio = sum(labels(:) ~= labels_old(:)) / num_voxels; fprintf(' 标签变化比例: %.4f\n', change_ratio); if change_ratio < convergence_thresh fprintf('3D算法在 %d 次迭代后收敛。\n', iter); break; end end % 清理标签 [labels, ~] = cleanup_labels_3d(labels); end function [clean_labels, map] = cleanup_labels_3d(labels) unique_vals = unique(labels(:)); map = containers.Map(unique_vals, 1:length(unique_vals)); clean_labels = zeros(size(labels)); for i = 1:numel(labels) clean_labels(i) = map(labels(i)); end end

3D版本注意事项：

• 数据顺序：Matlab中3D矩阵的维度顺序是(行, 列, 页)，通常对应(y, x, z)。在医学图像中，常表示为(height, width, depth)。代码中做了明确假设，使用时需与你的数据实际存储顺序保持一致。
• 计算开销：3D网格和距离计算内存消耗大。window_size不宜设置过大，否则window_mask计算会非常慢。对于大型3D数据，可能需要更优化的实现（如使用KD树搜索近邻）。
• 可视化：3D结果的可视化比2D复杂，通常需要查看切片或使用体绘制工具。

5. 完整示例：从加载图像到结果可视化

现在，我们将编写两个主脚本，分别演示2D和3D图像的分割全过程。

5.1 2D图像分割演示

创建一个名为main_2d_demo.m的脚本文件。

%% 基于局部质心的2D图像分割演示 clear; close all; clc; % 1. 加载或生成测试图像 % 示例1：使用Matlab内置图像 I = imread('cameraman.tif'); % 经典的摄影师图像 I = im2double(I); % 转换为double类型，范围[0,1] % 示例2：加载自己的图像 % I = im2double(imread('your_image.png')); % 2. 添加一些高斯噪声，模拟真实情况（可选） noise_level = 0.05; I_noisy = imnoise(I, 'gaussian', 0, noise_level^2); input_image = I_noisy; % 使用带噪声的图像作为输入 % 3. 显示原始图像 figure('Position', [100, 100, 1200, 400]); subplot(1,3,1); imshow(input_image); title('原始输入图像 (带噪声)'); colorbar; % 4. 设置算法参数 window_radius = 5; % 局部窗口半径。窗口大小为 (2*radius+1)^2。值越大，分割区域越大、越平滑。 max_iterations = 50; % 最大迭代次数 convergence_threshold = 0.001; % 收敛阈值。当标签变化比例小于此值时停止。 % 5. 调用分割函数 fprintf('开始2D局部质心分割...\n'); tic; [segmented_labels, final_centroids] = local_centroid_seg_2d(input_image, window_radius, max_iterations, convergence_threshold); elapsed_time = toc; fprintf('分割完成，耗时 %.2f 秒。\n', elapsed_time); % 6. 显示分割结果（标签图） subplot(1,3,2); imagesc(segmented_labels); axis image; title('分割标签图'); colormap(jet); % 使用jet色彩映射以区分不同标签 colorbar; % 7. 显示分割边界叠加图 % 计算标签图的边界 boundary_mask = boundarymask(segmented_labels); % 需要Image Processing Toolbox subplot(1,3,3); imshow(input_image); hold on; visboundaries(boundary_mask, 'Color', 'r', 'LineWidth', 1.5); % 红色边界叠加 title('分割边界叠加'); hold off; % 8. 输出一些统计信息 num_regions = max(segmented_labels(:)); fprintf('分割出 %d 个独立区域。\n', num_regions); fprintf('质心坐标示例（前5个）：\n'); disp(final_centroids(1:min(5, size(final_centroids,1)), :)); % 9. 保存结果（可选） % imwrite(uint8(255 * mat2gray(segmented_labels)), '2d_segmentation_result.png'); % save('2d_segmentation_data.mat', 'segmented_labels', 'final_centroids');

5.2 3D图像分割演示

创建一个名为main_3d_demo.m的脚本文件。由于3D数据较大，我们使用一个合成体积数据作为示例。

%% 基于局部质心的3D图像分割演示 clear; close all; clc; % 1. 生成合成3D体积数据（两个亮度不同的球体） [height, width, depth] = deal(64, 64, 64); % 创建一个64x64x64的小体积，便于快速演示 V = zeros(height, width, depth, 'double'); % 球体1参数 center1 = [height*0.3, width*0.3, depth*0.3]; radius1 = 15; intensity1 = 0.8; % 球体2参数 center2 = [height*0.7, width*0.7, depth*0.7]; radius2 = 12; intensity2 = 0.4; % 创建坐标网格 [X, Y, Z] = meshgrid(1:width, 1:height, 1:depth); % 生成球体1 dist1 = sqrt((Y - center1(1)).^2 + (X - center1(2)).^2 + (Z - center1(3)).^2); V(dist1 <= radius1) = intensity1; % 生成球体2 dist2 = sqrt((Y - center2(1)).^2 + (X - center2(2)).^2 + (Z - center2(3)).^2); V(dist2 <= radius2) = intensity2; % 添加一些高斯噪声 noise_level = 0.1; V = V + noise_level * randn(size(V)); V = max(0, min(1, V)); % 裁剪到[0,1]范围 fprintf('合成3D体积数据已创建。大小: %d x %d x %d\n', height, width, depth); % 2. 显示中间切片 figure; slice_idx = round(depth/2); subplot(1,2,1); imagesc(V(:,:,slice_idx)); axis image; colorbar; colormap(gray); title(sprintf('原始3D数据 (Z=%d切片)', slice_idx)); % 3. 设置算法参数（3D计算量更大，参数需更保守） window_radius = 3; % 3D窗口，半径3意味着7x7x7的立方体，计算量已不小 max_iterations = 30; convergence_threshold = 0.005; % 4. 调用3D分割函数 fprintf('开始3D局部质心分割...（这可能需要一些时间）\n'); tic; [segmented_labels_3d, final_centroids_3d] = local_centroid_seg_3d(V, window_radius, max_iterations, convergence_threshold); elapsed_time_3d = toc; fprintf('3D分割完成，耗时 %.2f 秒。\n', elapsed_time_3d); % 5. 显示分割结果的中间切片 subplot(1,2,2); imagesc(segmented_labels_3d(:,:,slice_idx)); axis image; colorbar; colormap(jet); title(sprintf('3D分割标签图 (Z=%d切片)', slice_idx)); % 6. 显示3D分割结果的等值面（可选，需要图形处理能力） % 提取一个主要区域（例如标签为2的区域）进行3D可视化 if license('test', 'Symbolic_Toolbox') % 简单检查3D绘图能力 try figure; label_to_visualize = 2; iso_value = 0.5; % 创建该标签的二进制体积 binary_volume = (segmented_labels_3d == label_to_visualize); % 使用isosurface进行渲染 patch(isosurface(binary_volume, iso_value), ... 'FaceColor', [0.8, 0.2, 0.2], 'EdgeColor', 'none', 'FaceAlpha', 0.6); daspect([1,1,1]); view(3); axis tight; camlight; lighting gouraud; title(sprintf('3D分割结果可视化 (标签 %d)', label_to_visualize)); xlabel('X'); ylabel('Y'); zlabel('Z'); grid on; catch fprintf('3D等值面渲染失败，可能由于数据量或图形支持问题。\n'); end end % 7. 输出统计信息 num_regions_3d = max(segmented_labels_3d(:)); fprintf('3D分割出 %d 个独立区域。\n', num_regions_3d); fprintf('3D质心坐标示例（前3个）：\n'); disp(final_centroids_3d(1:min(3, size(final_centroids_3d,1)), :)); % 8. 保存结果（可选） % save('3d_segmentation_results.mat', 'segmented_labels_3d', 'final_centroids_3d', 'V');

6. 运行结果与效果分析

运行main_2d_demo.m，你将看到类似下图的结果： （左侧）带噪声的原始图像 -> （中间）算法计算出的标签图，不同颜色代表不同区域 -> （右侧）分割边界（红色）叠加在原始图像上。

对于cameraman.tif这样的图像，算法能有效地将摄影师、背景和三角架分离，即使存在噪声。window_radius参数控制着分割的粒度。调大它，你会得到更少、更大的区域；调小它，则会得到更多、更细碎的分割。

运行main_3d_demo.m，你将在命令行看到迭代过程，并最终显示两个切面图和一个3D渲染图（如果系统支持）。算法成功地将两个亮度不同的3D球体分割开来。这证明了该方法在3D空间的有效性。

如何判断成功？

1. 视觉检查：分割边界是否大致贴合你感兴趣的物体边缘？区域内部是否均匀？
2. 收敛性：查看命令行输出的“标签变化比例”。一个健康的运行过程应该是该比例随着迭代迅速下降，并在若干次迭代后低于阈值，提示“算法已收敛”。
3. 区域数量：分割出的区域数量是否合理？对于简单图像，区域不应过多或过少。
4. 运行时间：对于2D图像（如512x512），算法应在数秒到数十秒内收敛。对于3D数据，时间会显著增加，这是迭代算法的固有特性。

7. 关键参数调优与常见问题排查

算法的表现很大程度上取决于参数设置。以下是详细的调优指南和问题排查表。

7.1 核心参数解析

7.2 常见问题、原因与解决方案

8. 工程实践与进阶建议

将算法用于实际项目时，以下几点能帮你走得更远：

1. 与经典方法对比：

   • K-means聚类：同样是无监督，但K-means基于全局灰度直方图，对空间信息不敏感，容易受噪声和灰度不均匀影响。局部质心法因其空间约束，通常对噪声和缓慢变化的背景更鲁棒。
   • 分水岭变换：对梯度敏感，极易导致过分割。通常需要结合标记点。局部质心法产生的结果通常比原始分水岭更紧凑。
   • Mean Shift：在概念上最接近。局部质心法可以看作是Mean Shift在规则图像网格上的一种特殊且高效的实现。

2. 处理复杂图像：

   • 纹理图像：纯灰度算法对纹理不敏感。可以考虑先用Gabor滤波器、LBP等方法提取纹理特征，将多特征图作为输入（需要修改权重计算方式）。
   • 彩色图像：将RGB图像转换到Lab或HSV颜色空间，对每个通道分别计算“质量”或使用欧氏距离融合多通道信息。
   • 医学图像（如MRI）：这类图像常存在强度不均匀性。局部质心法对此有一定抵抗力，但最好先进行偏置场校正。

3. 性能优化：

   • 多尺度策略：先在小尺度（下采样图像）上得到粗分割，再上采样到大尺度作为初始标签进行细化。这能大幅加速收敛并改善全局一致性。
   • 并行计算：算法的迭代过程天然可并行。可以使用Matlab的parfor循环来并行化“标签更新”步骤，特别是对于3D数据。
   • 使用编译语言重写核心循环：如果速度是瓶颈，可以考虑用C/C++编写MEX函数来计算距离和更新标签，在Matlab中调用。

4. 结果后处理与评估：

   • 区域合并：根据区域面积、平均灰度、形状等特征，将过分割的小区域进行合并。
   • 连通组件分析：使用bwconncomp或regionprops对二值化后的单个标签区域进行分析，过滤掉面积过小的噪声区域。
   • 无监督评估：在没有真实标注的情况下，可以使用内部指标评估分割质量，如区域内部的灰度一致性（方差小为好）和区域之间的灰度对比度（差异大为好）。

基于局部质心的无监督分割是一个强大而直观的工具箱。它完美诠释了如何将简单的物理概念转化为有效的图像分析算法。通过本文，你不仅获得了一套即拿即用的Matlab代码，更重要的是理解了其内在机理和调参逻辑。下次当你面对没有标签的数据时，不妨先尝试这个“物理驱动”的方法，它可能会给你带来一个干净、可靠的起点。建议将本文代码收藏，作为你图像处理武器库中的一件灵活备选方案。