Kmeans算法和FCM算法(含MATLAB代码和具体实验)

Kmeans算法和FCM算法

算法原理

    以上概念都是基于图像方面给出的，因为之前做过基于K均值改进算法的图像方面的研究，所以我就直接套用文章中的内容了，算法的原理就是根据目标范围中考察点之间的距离来分类的，而像素在计算机中就是一个个坐标，这也是为什么算法能应用于图像的原因。

例子

    这个图像来源于本站某位大佬的算法运行截图（一时间没找到出处，望海涵），左图是原始点位分布，右图是以4为聚类数的聚类效果图（图中黑色的正方形就是聚类中心）。

代码

Kmeans算法

function [idx,C] = my_kmeans(X,k)
% 输入:X为数据,k为分类数
% 输出:idx为类的索引,C为聚类好的聚类中心
% matlab自带kmeans函数，调用格式为"[idx,C] = kmeans(X,k)",
% 因为是用k-means++算法进行簇中心初始化，所以比我写的算法性能要好
%% 原理推导Kmeans聚类算法
[m,n]=size(X);
C=X(randperm(m,k),:);%从m个点中随机选择cluster_num个点作为初始聚类中心点
iter_max=100;%最大次数
d_C = 1e-5; %中心变化阈值
iter=0;
distance1 = zeros(k,m);
cluster_new = zeros(k,n);
while(iter<iter_max)
    % distance1存储每个点到各聚类中心的欧氏距离
    for i=1:k
        distance=(X-C(i,:)).^2;% 
        distance1(i,:)=sqrt(sum(distance')); % sqrt(x^2 + y^2),distance1是50行i列的
    end
    [~,idx]=min(distance1);%index_cluster取值范围1~cluster_num,~是最小元素的值,用不到,故省略
    % cluster_new存储新的聚类中心
    for j=1:k
        cluster_new(j,:)=mean(X(idx==j,:));
    end
    %如果新的聚类中心和上一轮的聚类中心距离和大于therad_lim，更新聚类中心，否则算法结束
    if (sum(sqrt(sum(((cluster_new-C).^2)')))>d_C)
        C=cluster_new;
    else
        break; 
    end
    iter = iter+1;
end
fprintf('迭代次数:%d\n',iter)
idx = idx';%为了和MATLAB自带算法输出维度一样
end

说明：MATLAB翻译过来就是“矩阵实验室”，该软件对矩阵类数据处理起来很快，所以在学习使用MATLAB时，尽量要将算法写成矩阵运算的形式，而不是用for循环一个一个的去算元素，虽然都可以实现算法的功能，但是一旦数据量上去了，不同的代码形式会产生十分明显的运算时间差异。

Kmeans++算法

    在Kmeans算法代码注释中我提到了MATLAB自带Kmeans算法，其实自带的是Kmeans++算法，调用格式为

[idx,C] = kmeans(X,k)

可以用

help kmeans

了解算法
    Kmeans++相对于Kmeans就是改进了初始化中心的方式，让中心在初始化的时候不那么随机，而是要相隔一段距离，这样做是为了避免初始化的时候中心隔的太近，造成最终聚类效果差。

FCM算法

function [center,U] = my_fcm(data,Nc)
% 输入:data为数据,Nc为分类数
% 输出:centers为聚类好的聚类中心,U为隶属度矩阵
% matlab自带fcm函数，调用格式为"[centers,U] = fcm(data,Nc)",
%%
max_iter = 100;%迭代次数
m = 2;%指数
n = size(data,1);%样本个数
%--初始化隶属度u，条件是每一列和为1
U = rand(Nc,n);
col_sum = sum(U);
U = U./col_sum(ones(Nc,1),:);
% d_U = 1e-5;         % 隶属度变化阈值
d_C = 1e-5;         % 类中心变化阈值
%% 
%根据算法编写的代码，运行效率太低了，得改用矩阵运算的格式
% for i = 1:iter
%     %更新centers
%     for j = 1:Nc
%         u_ij_m = U(j,:).^m;
%         sum_u_ij = sum(u_ij_m); 
%         centers(j,:) = u_ij_m*data./sum_u_ij;
%     end
%     %计算目标函数J
%     temp1 = zeros(Nc,n);
%     for j = 1:Nc
%         for k = 1:n
%             temp1(j,k) = U(j,k)^m*(norm(data(k,:)-centers(j,:)))^2;
%         end
%     end
%     J(i) = sum(sum(temp1));
%     %更新U
%     for j = 1:Nc
%         for k = 1:n
%             sum1 = 0;
%             for j1 = 1:Nc
%                 temp = (norm(data(k,:)-centers(j,:))/norm(data(k,:)-centers(j1,:))).^(2/(m-1));
%                 sum1 = sum1 + temp;
%             end
%             U(j,k) = 1./sum1;
%         end
%     end
%     i
% end

iter = 0;
mf = U.^m;
center = mf*data./((ones(size(data, 2), 1)*sum(mf'))');  
while iter < max_iter

    dist = distfcm(center, data);        % 计算距离矩阵
    tmp = dist.^(-2);
    U_new = tmp./(ones(Nc, 1)*sum(tmp));	% 新的隶属度矩阵
    mf = U_new.^2;                              % 隶属度矩阵进行指数运算
    center_new = mf*data./((ones(size(data, 2), 1)*sum(mf'))'); % 新聚类中心

%     err_U = U - U_new;
    err_C = center - center_new;
    
%     if(max(abs(err_C(:))) < d_C || max(abs(err_U(:))) < d_U)
    if(max(abs(err_C(:))) < d_C )
        break;
    end
    U = U_new;
    center = center_new;
    iter = iter+1;
end
fprintf('迭代次数:%d\n',iter)
end

    FCM算法相对于Kmeans算法，加入了模糊的概念，对目标函数也进行相应的修改，通过伟大的拉格朗日乘法求解，得到隶属度函数的表达式，具体的求解过程本站也有大佬推理过。
    MATLAB自带FCM算法，调用格式为

[centers,U] = fcm(data,Nc)

实验

    以下两张图片都是含病变组织（脑瘤）的脑部MRI图，分别命名为“001.png”，“002.png”，放在mian函数的同一个文件夹中

实验代码

运用Kmeans算法对001图进行聚类分割

clear;clc;close all;
%% 图像导入和预处理
data_original = imread('001.png');
figure; imshow(data_original);
title('原始图像')
if (size(data_original,3)>2)
data_original = rgb2gray(data_original);
end
%将图像转换为双精度值
img = im2double(data_original);
%通过应用中值滤波器去除噪声
img=medfilt2(img);
figure; imshow(img);
title('去噪后的图像')

%% 用Kmeans算法分割
%将图像转换为线性形状
data = reshape(img,[],1);

Nc = 4;%聚类数

%使用Kmeans方法
 [idx,center] = my_kmeans(data,Nc);% [idx,center] = kmeans(data,Nc);自编函数和自带函数都能用
%聚类后的图
im_res = reshape(center(idx),size(img));
figure;imshow(im_res);
img_res = reshape(idx,size(img));
figure;imagesc(im_res);

%分割后的图，这里是4类，所以是四张图
figure
subplot(2,2,1),imshow(img_res==1,[]);
subplot(2,2,2),imshow(img_res==2,[]);
subplot(2,2,3),imshow(img_res==3,[]);
subplot(2,2,4),imshow(img_res==4,[]);

运用FCM算法对002图进行聚类分割

clear;clc;close all;
%% 图像导入和预处理
data_original = imread('002.png');
figure; imshow(data_original);
title('原始图像')
if (size(data_original,3)>2)
data_original = rgb2gray(data_original);
end
%将图像转换为双精度值
img = im2double(data_original);
%通过应用中值滤波器去除噪声
img=medfilt2(img);
figure; imshow(img);
title('去噪后的图像')

%% 用FCM算法分割
%将图像转换为线性形状
data = reshape(img,[],1);

Nc = 4;%聚类数

%使用FCM或AFKM
[center,U] = my_fcm(data,Nc);%[center,U] = fcm(data,Nc);自编函数和自带函数都能用
[maxU,idx] = max(U);
%聚类后的图
im_res = reshape(center(idx),size(img));
figure;imshow(im_res);
title('聚类后的图')
%聚类后的图(一种颜色代表一类)
img_res = reshape(idx,size(img)); 
figure;imagesc(img_res)
title('聚类后的彩图')

%分割后的图，这里是4类，所以是四张图
figure
subplot(2,2,1),imshow(img_res==1,[]);
subplot(2,2,2),imshow(img_res==2,[]);
subplot(2,2,3),imshow(img_res==3,[]);
subplot(2,2,4),imshow(img_res==4,[]);

实验结果

    对于001图

    对于002图

对002图展开说明：


下面的步骤可以提取病变组织部分，但是不属于本文章所提供的内容