图像分割之k-means聚类分割

k-means聚类原理

k-means聚类是一种无监督学习的聚类算法，它的目的是将数据集中的样本划分成若干个类别，使得同一类别内的样本相似度高，而不同类别之间的样本相似度低。其目标是将数据点划分为k个类簇，找到每个簇的中心并使其度量最小化。
此外，K-means算法有一些优点，如简单易懂、计算速度快，适合处理大数据集。但它也有一些缺点，比如需要预先指定簇的数量K，对初始聚类中心的选择敏感，可能收敛到局部最优而非全局最优解，对于非凸数据集和噪声数据敏感。
在实际应用中，K-means算法广泛应用于市场细分、图像分割、社交网络分析等领域。本文主要将其应用于图像分割领域。

k-means基本思路

通过计算相似度（默认欧氏距离），将相似度大的样本聚集到同一个类别，k表示聚成k个类别，means表示每个类别的聚类中心点是通过簇中所有样本点的均值得到。

聚类效果评估

聚类的效果可采用误差平方和SSE进行评价，SSE值越小，表示数据点越接近它们的中心点，聚类效果越好。

k-means聚类算法流程

下面是K-Means聚类算法的分析流程，步骤如下：

第一步，确定K值，即将数据集聚集成K个类簇或小组。
第二步，从数据集中随机选择K个数据点作为质心（Centroid）或数据中心。
第三步，分别计算每个点到每个质心之间的距离，并将每个点划分到离最近质心的小组，跟定了那个质心。
第四步，当每个质心都聚集了一些点后，重新定义算法选出新的质心。
第五步，比较新的质心和老的质心，如果新质心和老质心之间的距离小于某一个阈值，则表示重新计算的质心位置变化不大，收敛稳定，则认为聚类已经达到了期望的结果，算法终止。
第六步，如果新的质心和老的质心变化很大，即距离大于阈值，则继续迭代执行第三步到第五步，直到算法终止。

对于图像分割任务，需要在第一步之前将图像reshape为m*n行、p列的矩阵，这样每一行就是一个p维数据点。通过以上聚类算法，得到最终的聚类结果，最后再将结果reshape为原始图像的大小即可。

MATLAB实现

MATLAB自带函数

kmeans函数的说明如下：

idx = kmeans(X,k) 执行 k 均值聚类，以将 n×p 数据矩阵 X 的观测值划分为 k个聚类，并返回包含每个观测值的簇索引的 n×1 向量 (idx)。X 的行对应于点，列对应于变量。

根据上面的kmeans函数的说明，我们需要将图像矩阵reshape一下，然后再调用该函数，最后再reshape回来，下面是测试的主要代码：

[m,n,p]=size(img)              %m,n为所求，p=3为通道数
% 将图像进行RGB——3通道分解
% org(:, :, 1)......分别代表rgb通道
A = reshape(img(:, :, 1), m*n, 1);   
B = reshape(img(:, :, 2), m*n, 1);
C = reshape(img(:, :, 3), m*n, 1);
data = [A B C];  % r g b分量组成样本的特征，每个样本有三个属性值，共width*height个样本

% kmeans第一个参数: N*P的数据矩阵，N为数据个数，P为单个数据维度
res = kmeans(double(data), k);        
result = reshape(res, m, n);  %转化为图片形式
1
2
3
4
5
6
7
8
9
10
11

MATLAB实现

根据上述算法流程，使用matlab编写如下代码，其中设置强制截断的迭代次数，并添加了迭代过程展示。
my_kmeans函数：

function [C, label, J] = my_kmeans(data, k, is_gpu, is_show) 
%% 语法 ：
% [C, label, J] = my_kmeans(data, k, is_gpu, is_show) 
%% 说明：
% data  为输入图像
% k  为聚类的k个数据中心
% is_gpu  是否启用gpu加速
% is_show  参数决定是否显示聚类过程
% ----------------------------------------
% return label  为返回的聚类结果
% return C  聚类中心
% return J  每次迭代的误差
%%
if nargin >= 3
    is_gpu = 1;
else
    is_gpu = 0;
end
if nargin == 4
    is_show = 1;
else
    is_show = 0;
end
%%
% data可以为一维数据
m = size(data, 1);
n = size(data, 2);
p = size(data, 3);
% [m, n, p] = size(I); %图片的大小m*n，p代表RGB三层
X = reshape(double(data), m*n, p);
if is_gpu
    X = gpuArray(X);
end
rng('default');
C = X(randperm(m*n, k), :);%随机选k个聚类中心

J = [];
J_prev = inf; 
tol = 1e-1; %容忍度tol，inf为无穷大
iter = 0; 
rand_c = rand(k,3)*255;
rand_c = uint8(rand_c);
while true
    iter = iter + 1;
    dist = sum(X.^2, 2)*ones(1, k) + (sum(C.^2, 2)*ones(1, m*n))' - 2*X*C'; %计算各个点到K个聚类中心的距离
    [~,label] = min(dist, [], 2) ;  %label记录最小值的行数
    for i = 1:k
       C(i, :) = mean(X(label == i , :)); %取新的k个聚类中心
    end
    J_cur = sum(sum((X - C(label, :)).^2, 2)); %距离之和
    J = [J, J_cur];
    fprintf('#iteration: %03d, objective norm diff: %f\n', iter, norm(J_cur-J_prev, 'fro'));
    if norm(J_cur-J_prev, 'fro') < tol   % A和A‘的积的对角线和的平方根，即sqrt(sum(diag(A'*A)))，本次与上次距离之差
        break;
    end
    if (iter==100)  %设置强制截断的迭代次数
        break;
    end
    J_prev = J_cur;
    if is_show
        % 录制gif
        rand_c = uint8(C);
        temp = uint8(label);
        temp = rand_c(temp,:);
        sg=reshape(temp,m,n, []);
        imshow(mat2gray(sg));
        F=getframe(gcf);
        data=frame2im(F);
        [data,map]=rgb2ind(data,256);
        if iter == 1
            imwrite(data,map,'test.gif','gif','Loopcount',inf,'DelayTime',0.2);
        else
            imwrite(data,map,'test.gif','gif','WriteMode','append','DelayTime',0.2);
        end
    end
end
% C = rand_c;
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78

测试代码如下，如果是单纯是为了处理图像，可以将“聚类”这小段的代码重新封装成一个function。

clear;
close all;
clc;
tic
%% 读取数据
img=imread('city.png');
% img = rgb2gray(img);
%% 聚类
m = size(img, 1);
n = size(img, 2);
p = size(img, 3);
k = 4;
[C, label, J2] = my_kmeans(img, k, 1, 1);
% 使用gather函数将gpuarray转换为普通数组
C = gather(C);
C = uint8(C);
% randn('seed', sum(100*clock))
% sprintf('seed=%d', sum(100*clock))
% rand_c = rand(k,3)*255;
% C = uint8(rand_c);
label = uint8(label);
temp = C(label,:);
dst = reshape(temp, m, n, []);%转换成图片矩阵的格式
toc
%% 显示结果
figure
subplot(221), imshow(img,[]);
subplot(222), imshow(dst,[]);
subplot(212), plot(J2), xlabel('iters'), ylabel('norm diff')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

测试结果

测试的原始图片如下：

聚类过程如下，是不是有种延时摄影的感觉。

测试结果及相应的迭代过程如下所示：

MATLAB自带函数不同K值下的效果：

自己实现的不同K值下的效果：

MATLAB自带函数的效果和我实现的效果有所区别，主要是由于MATLAB自带函数采用了kmeans++进行了优化，在此就不赘述了，感兴趣的朋友可以去深入了解。

参考文献

[1] 机器学习实战——K-means聚类图像分割
[2] 基于K均值聚类算法的图像分割（Matlab）