简介：这份资源是Normalized Cuts（规范化切割）算法的Matlab实现代码包，面向计算机视觉、图像处理方向的学习者与研究者，帮助理解并复现这一经典图论分割方法。包内共77个文件，以50个m脚本、12个cpp源码、5个mexw64编译文件为主，另有少量jpg示例图、txt说明与db数据文件，压缩包约103KB，涵盖相似度矩阵构建、归一化拉普拉斯矩阵计算、特征分解、特征向量选择与K-means聚类等完整流程，并附有图像分割与聚类演示脚本。已有344人学习下载。读者可借助这些代码深入理解Normalized Cuts从图模型构建到谱聚类求解的实现细节，对照示例图像验证分割效果，并在此基础上扩展到视频分割、社区发现等场景，是入门谱聚类与图像分割的实用参考。

1. 从一次图像分割翻车说起：这个 Matlab 实现到底能干什么

去年帮一个做医学影像的朋友处理一批细胞切片图，他丢过来一句“把前景背景分开就行”，我心想阈值分割三行代码搞定。结果光照不均加上细胞粘连，阈值法直接翻车，边缘碎得像被狗啃过。后来翻出这份 Normalized Cuts 的 Matlab 实现，才把粘连区域干净地切开。Normalized Cuts（简称 Ncut）是谱聚类在图像分割里的经典落地，核心思路是把图像看成一张带权图，像素是节点，像素间的相似度是边权，然后通过求解广义特征值问题把图切成若干块，让“块内相似度高、块间相似度低”。这份资源就是它在 Matlab 环境下的完整实现，适配 Matlab2014a 及以上版本，适合做图像分割、区域提取、聚类实验的开发者，尤其是需要快速验证算法效果又不想从零推导公式的人。它不依赖深度学习框架，纯矩阵运算，拿过来改改相似度函数就能跑自己的数据。

2. 算法骨架与工程结构：先看清 Ncut 是怎么被拆成代码的

2.1 从图论到矩阵：Ncut 的数学骨架

Ncut 的目标函数长这样：给定图 G=(V,E)，把节点分成 A、B 两部分，最小化

Ncut(A,B) = cut(A,B)/assoc(A,V) + cut(A,B)/assoc(B,V)

其中 cut(A,B) 是 A 和 B 之间所有边的权重和，assoc(A,V) 是 A 中节点与全图所有边的权重和。这个归一化项是关键，它避免了最小割偏向切掉孤立小区域。直接优化是 NP 难的，所以工程上把它松弛成连续问题：求解 (D-W)x = λDx 的广义特征值问题，W 是相似度矩阵，D 是度矩阵。取第二小特征值对应的特征向量，用它做阈值分割或 k-means 聚类，就得到最终分割结果。

这份 Matlab 实现就是围绕这个流程展开的。它没有用任何工具箱里的谱聚类函数，而是手写相似度矩阵构建、稀疏矩阵求解和特征向量后处理，方便你逐行看懂每一步在干什么。

2.2 目录结构与模块职责

拿到源码包，目录大致是这样组织的（不同版本可能略有差异，但核心文件跑不掉）：

Ncut/ ├── main.m % 入口脚本，读图、调参、显示结果 ├── computeW.m % 构建相似度矩阵 W ├── computeD.m % 由 W 计算度矩阵 D ├── ncut.m % 核心：求解广义特征值，返回分割向量 ├── discretize.m % 把连续特征向量离散化成标签 ├── imgsegment.m % 封装完整流程，输入图像输出分割图 ├── utils/ % 辅助函数 │ ├── imread_color.m % 统一读图接口 │ └── showmask.m % 叠加显示分割掩膜 └── data/ % 示例图片

main.m是给你直接跑的，里面通常预设了一张示例图。computeW.m和computeD.m负责图构建，ncut.m是计算最重的部分，discretize.m决定最终标签怎么定。这种拆分方式很“教学向”，每个文件对应一个数学步骤，改起来定位快。

2.3 相似度矩阵构建：W 的三种常见写法

相似度矩阵 W 直接决定分割质量。这份实现里computeW.m一般提供两种模式：全连接和稀疏连接。全连接就是每个像素和所有像素算相似度，公式是

W(i,j) = exp(-||F(i)-F(j)||^2 / sigma_F) * exp(-||X(i)-X(j)||^2 / sigma_X)

F 是像素特征（灰度或颜色），X 是空间坐标。两个 sigma 分别控制特征相似度和空间距离的敏感度。全连接矩阵是 N×N 稠密矩阵，N 是像素数，一张 200×200 的图就是 40000×40000，内存直接爆炸。所以实际代码里通常用稀疏化：只保留每个像素周围 r 邻域内的连接，或者用 kd-tree 找 k 近邻。

% computeW.m 片段：稀疏邻域相似度构建 function W = computeW(img, sigmaF, sigmaX, radius) [h, w, ~] = size(img); N = h * w; % 把图像拉成 N x 3 的特征矩阵（RGB） feat = reshape(img, N, 3); % 生成像素坐标 [xx, yy] = meshgrid(1:w, 1:h); coord = [xx(:), yy(:)]; % 用 kd-tree 找每个点的邻域，避免全连接 % 常见做法是 rangesearch 或 knnsearch idx = rangesearch(coord, coord, radius); % 逐点计算权重，组装稀疏矩阵 I = []; J = []; S = []; for i = 1:N nb = idx{i}; for j = nb if i == j, continue; end df = sum((feat(i,:) - feat(j,:)).^2); dx = sum((coord(i,:) - coord(j,:)).^2); w = exp(-df / sigmaF) * exp(-dx / sigmaX); I(end+1) = i; J(end+1) = j; S(end+1) = w; end end W = sparse(I, J, S, N, N); W = (W + W') / 2; % 保证对称 end

这段代码里sigmaF和sigmaX是调参重点。sigmaF越大，颜色差异被容忍得越多，分割区域更倾向于按空间连通性走；sigmaX越大，空间距离影响越弱，容易出现远距离同色区域被合并。radius控制邻域半径，一般取 5 到 15 像素。注意最后做了对称化，因为逐点计算可能因为浮点误差导致 W 不对称，而后续特征值求解要求对称矩阵。

2.4 广义特征值求解：ncut.m 里的计算核心

ncut.m是整个项目最“重”的文件。它要解的是

(D - W) x = λ D x

Matlab 里用eigs求解前 k 个最小特征值。因为 D 是对角矩阵，可以转成标准特征值问题：令 y = D^(1/2) x，则

D^(-1/2) (D - W) D^(-1/2) y = λ y

这样就能用eigs直接算对称矩阵的特征值。

% ncut.m 片段：求解第二小特征向量 function [eigvec, eigval] = ncut(W, D, k) N = size(W, 1); % 构造归一化拉普拉斯矩阵 Dinv_sqrt = spdiags(1./sqrt(diag(D)), 0, N, N); L = Dinv_sqrt * (D - W) * Dinv_sqrt; L = (L + L') / 2; % 强制对称，避免数值误差 % 求前 k 个最小特征值对应的特征向量 % 注意：最小特征值恒为 0，对应平凡解，要跳过 [V, E] = eigs(L, k + 1, 'sm'); eigval = diag(E); eigvec = V(:, 2:k+1); % 丢掉第一个平凡解 % 把特征向量还原回原空间 eigvec = Dinv_sqrt * eigvec; end

这里有几个坑点。第一，eigs的'sm'选项在小矩阵上可能不稳定，如果报错可以换成'sa'或者直接eig(full(L))，但后者对大规模矩阵不现实。第二，求出来的特征向量需要归一化，常见做法是按行做 L2 归一化，让每个像素的特征向量长度为 1，这样后续聚类对尺度不敏感。第三，k是你要分割的区域数，二分割取 k=1，多分割取 k>1。代码里eigs多求一个特征值就是为了丢掉那个恒为 0 的平凡解。

2.5 离散化：从连续特征向量到最终标签

特征向量是连续的，不能直接当标签用。discretize.m负责把它变成离散区域编号。二分割最简单：取第二小特征向量，找一个阈值（常用 0 或者中位数）把像素分成两类。多分割则要对前 k 个特征向量做 k-means 聚类。

% discretize.m 片段：二分割与多分割 function labels = discretize(eigvec, k) if k == 1 % 二分割：用零阈值或中位数 v = eigvec(:, 1); thresh = median(v); % 比 0 更鲁棒 labels = double(v > thresh) + 1; else % 多分割：对特征向量矩阵做 k-means % 常见做法是先按行归一化 X = eigvec ./ repmat(sqrt(sum(eigvec.^2, 2)), 1, k); labels = kmeans(X, k, 'Replicates', 5, 'MaxIter', 200); end end

二分割里用中位数而不是 0 做阈值，是因为特征向量的分布不一定以 0 为中心，中位数能保证两类像素数量大致均衡。多分割里kmeans的Replicates参数设成 5 是为了避免陷入局部最优，代价是计算时间翻几倍。如果数据量大，可以降到 2 或 3。

3. 跑起来：从零到一张分割图的完整操作

3.1 环境准备与版本适配

这份实现针对 Matlab2014a 及以上版本。2014a 是个分水岭，之前的版本eigs对稀疏矩阵支持不够好，kmeans的参数名也有差异。如果你用的是 2016b 之后的版本，基本可以直接跑，但要注意两点：一是rangesearch在旧版本里叫rangesearch，新版本没变，但knnsearch的返回参数格式在 2019a 之后有调整；二是kmeans的Replicates参数在 2014a 里是支持的，但默认距离度量是sqeuclidean，新版本默认也是这个，一般不用改。

没有额外工具箱依赖，Image Processing Toolbox 只在读图和显示时用到，核心算法纯基础矩阵运算。如果你连 IPT 都不想装，把imread换成imread原生函数、显示部分用imagesc替代即可。

3.2 最小可运行示例：跑通 main.m

打开main.m，通常长这样：

% main.m：最小可运行入口 clear; clc; close all; % 1. 读图 img = imread('data/sample.jpg'); img = imresize(img, [128, 128]); % 降采样，控制矩阵规模 figure; imshow(img); title('原图'); % 2. 参数设置 sigmaF = 0.1; % 特征相似度尺度 sigmaX = 5; % 空间距离尺度 radius = 8; % 邻域半径 k = 2; % 分割区域数 % 3. 构建 W 和 D W = computeW(double(img)/255, sigmaF, sigmaX, radius); D = computeD(W); % 4. 求解 Ncut [eigvec, eigval] = ncut(W, D, k-1); % 5. 离散化 labels = discretize(eigvec, k-1); % 6. 显示结果 mask = reshape(labels, size(img,1), size(img,2)); figure; imshow(mask, []); title('分割结果');

跑之前先把data/sample.jpg换成你自己的图，或者把imread路径改成绝对路径。imresize到 128×128 是为了控制计算量，N=16384，稀疏矩阵还能扛住。如果你直接上 512×512，N=262144，eigs会跑很久，内存也吃紧。

3.3 参数怎么调：sigmaF、sigmaX、radius 的联动

这三个参数没有理论最优值，只能按图调。我一般按这个顺序试：

调参时先固定 radius=8，sigmaX=5，然后只动 sigmaF。如果分割结果把该分开的区域粘在一起，说明 sigmaF 太大，往下调；如果同一区域被切碎，说明 sigmaF 太小，往上调。sigmaX 一般不用大动，除非你的图像分辨率特别高或特别低。

3.4 结果可视化与掩膜叠加

showmask.m通常做一件事：把分割标签转成彩色掩膜，叠加到原图上。核心代码就几行：

% showmask.m：叠加显示 function showmask(img, labels) % 生成随机颜色表 nLabels = max(labels); colors = rand(nLabels, 3); % 把标签映射成 RGB maskRGB = reshape(colors(labels, :), size(img)); % 叠加：原图 0.6 + 掩膜 0.4 overlay = 0.6 * double(img)/255 + 0.4 * maskRGB; figure; imshow(overlay); title('分割叠加'); end

这样能直观看到分割边界落在哪里。如果边界和实际物体边缘对不齐，回去调 sigmaF 和 radius。

4. 避坑与排查：跑 Ncut 最容易翻车的五个地方

4.1 现象：eigs 报错“矩阵太大”或直接卡死

原因：全连接 W 是 N×N 稠密矩阵，N 超过 5000 就基本不可行。很多人第一次跑忘了稀疏化，或者 radius 设得太大导致稀疏矩阵退化成稠密。

解决：确认computeW.m里用的是sparse而不是zeros。把 radius 降到 5~8，或者先用imresize把图像降到 128×128 以下。如果还卡，改用eigs(L, k+1, 'sm', 'Tolerance', 1e-6, 'MaxIter', 300)限制迭代次数。

4.2 现象：分割结果全黑或全白，没有分界

原因：特征向量阈值选错了。二分割时如果直接用 0 做阈值，而特征向量整体偏移，就会全分到一类。

解决：把阈值改成median(v)或mean(v)。更稳的做法是对特征向量做归一化后再取中位数。如果还是不行，检查ncut.m里是不是忘了丢掉第一个平凡解，导致拿到的特征向量恒为常数。

4.3 现象：kmeans 报错“X 包含 NaN 或 Inf”

原因：特征向量里有 NaN，通常是因为 D 的对角线有 0，求逆时除了零。D 的对角线是 W 的行和，如果某个像素和所有像素的相似度都下溢到 0，行和就是 0。

解决：在computeD.m里加一个极小值保护：D = spdiags(max(sum(W,2), 1e-10), 0, N, N)。或者在ncut.m里对 Dinv_sqrt 做同样处理。这个坑在图像有大片纯黑或纯白区域时特别常见。

4.4 现象：分割边界呈块状锯齿，不平滑

原因：Ncut 本身是像素级分割，没有做任何平滑或后处理。加上降采样后，边界自然粗糙。

解决：如果只是可视化需求，可以在discretize之后加一个中值滤波或形态学闭运算。如果是精度需求，别降采样，直接在全分辨率上跑，但要做好等几分钟的准备。另一个办法是增大 radius，让邻域连接更密，边界会更连续，代价是计算量上升。

4.5 现象：换一张图就完全不能用，参数全要重调

原因：sigmaF 和 sigmaX 是绝对尺度，对图像亮度和分辨率敏感。一张图调好的参数，换一张图就失效。

解决：把 sigmaF 改成自适应：先算全图特征差值的标准差，令 sigmaF = 0.1 * std(df)。sigmaX 同理，用坐标差值的标准差做基准。这样参数变成相对值，跨图迁移性好很多。这个改动我一般直接写在computeW.m开头，省得每次手动调。

5. 二次开发：把 Ncut 改成你自己的分割工具

5.1 替换相似度函数：从 RGB 到多特征融合

默认的 W 只用 RGB 和坐标。如果你做医学影像，可以把灰度、纹理、梯度方向都拼进特征向量。改法很简单：在computeW.m里把feat从 N×3 扩成 N×d，然后调整 sigmaF 的维度。比如加一个 LBP 纹理特征：

% 扩展特征：RGB + LBP gray = rgb2gray(img); lbp = extractLBPFeatures(gray); % 需要 Computer Vision Toolbox % 或者手写简化版 LBP feat = [reshape(img, N, 3), reshape(lbp, N, 1)]; % sigmaF 要相应调大，因为维度变了

注意特征维度增加后，欧氏距离的尺度会变，sigmaF 要重新标定。我一般先用pdist算一下特征距离的分布，取中位数的 0.1 倍作为 sigmaF 初值。

5.2 加速策略：从 eigs 到 Lobpcg 或子采样

eigs对大规模稀疏矩阵还是慢。两个加速方向：一是用 LOBPCG 算法，Matlab 没有内置，但可以自己实现或者找第三方；二是子采样，先在小图上算特征向量，再用插值或最近邻把标签映射回全分辨率。第二种更实用，代码改动也小：

% 子采样加速：先降采样求标签，再上采样 img_small = imresize(img, 0.25); labels_small = imgsegment(img_small, sigmaF, sigmaX, radius, k); labels_full = imresize(labels_small, size(img(:,:,1)), 'nearest');

这样计算量降到 1/16，边界会粗一点，但大多数场景够用。如果边界精度要求高，可以在全分辨率上只对边界附近的像素重新计算。

5.3 多区域分割的标签合并技巧

k>2 时，kmeans 出来的标签顺序是随机的，同一张图跑两次可能标签编号对不上。如果你要做批量处理，需要固定标签顺序。常见做法是按区域面积排序，或者按区域质心的空间位置排序。我一般用面积降序：

% 固定标签顺序：按面积从大到小重排 areas = histcounts(labels, 1:k+1); [~, order] = sort(areas, 'descend'); map = zeros(1, k); map(order) = 1:k; labels = map(labels);

这样每次跑出来的标签 1 都是最大区域，方便后续统计和对比。

5.4 一个我踩过的坑：稀疏矩阵对称化的代价

早期我为了省事，在computeW.m里只算上三角，然后W = W + W'。结果发现分割结果左右不对称。排查半天才意识到，sparse组装时如果 I、J 有重复索引，Matlab 会把对应值相加而不是覆盖。上三角加转置后，对角线元素被加了两次。后来改成先组装完整矩阵再对称化，或者组装时跳过对角线。从那以后我每次写完稀疏矩阵都强制走一遍issymmetric(W)检查，不通过就不往下跑。希望这个习惯也能帮你省掉几个小时的排查时间。

本文还有配套的精品资源，点击获取