[1] Chuang K S, Tzeng H L, Chen S, et al. Fuzzy c-means clustering with spatial information for image segmentation[J]. computerized medical imaging and graphics, 2006, 30(1): 9-15.
概述模糊c均值聚类（FCM）算法用于图像分割时，由于缺少图像的空间信息，因此对噪声的鲁棒性较差，为此国内外学者提出了许多FCM变体算法。其中主流的方法是通过构造新的目标函数将空间信息添加到原始FCM中，从目前来看，构造目标函数的方法大致分为两类：一是通过滤波算法预处理，接着在处理后的图像上利用数学知识构造新的目标函数，达到提高鲁棒性的目的；二是利用局部信息或正则化运算改变隶属度矩阵，进而达到提高鲁棒性的目的。
本文主要关注的第二种方法，即通过局部运算，改变每次迭代的隶属度值，从而将图像的空间信息引入FCM目标函数之中，提高了算法的鲁棒性。题主根据论文中的公式，复现了这篇论文的代码，从效果来看，该算法的噪声平滑能力较为有限，但
V
PC
{{text{V}}_{text{PC}}}
VPC和
V
PE
{{text{V}}_{text{PE}}}
VPE的值却比原始论文的高。
- 参数设置：包括聚类数目 c l u s t e r _ n u m cluster_{num} cluster_num、模糊因子 m m m、最大迭代次数 m a x _ i t e r max_{iter} max_iter、窗口尺寸 k k k；
- 图像灰度转换及其归一化处理；
- 添加噪声（可无）；
- 初始化隶属度矩阵 u i j {{u}_{ij}} uij；
- 更新聚类中心 c i c_{i} ci，更新公式为 c i = ∑ j = 1 N u i j m x j ∑ j = 1 N u i j m {{c}_{i}}=frac{sumlimits_{j=1}^{N}{u_{ij}^{m}{{x}_{j}}}}{sumlimits_{j=1}^{N}{u_{ij}^{m}}} ci=j=1∑Nuijmj=1∑Nuijmxj；
- 更新隶属度矩阵 u i j u_{ij} uij，更新公式为 u i j = 1 ∑ k = 1 c ( ∥ x j − c i ∥ ∥ x j − c k ∥ ) 2 / ( m − 1 ) {{u}_{ij}}=frac{1}{sumlimits_{k=1}^{c}{{{left( frac{left| {{x}_{j}}-{{c}_{i}} right|}{left| {{x}_{j}}-{{c}_{k}} right|} right)}^{{2}/{left( m-1 right)};}}}} uij=k=1∑c(∥xj−ck∥∥xj−ci∥)2/(m−1)1；
- 计算空间函数 h i j h_{ij} hij，计算公式为 h i j = ∑ k ∈ N ( x j ) u i k {{h}_{ij}}=sumlimits_{kin Nleft( {{x}_{j}} right)}{{{u}_{ik}}} hij=k∈N(xj)∑uik；
- 计算新的隶属度矩阵 u ′ i j {{{u}'}_{ij}} u′ij，计算公式为 u ′ i j = u i j p h i j q ∑ k = 1 c u k j p h k j q {{{u}'}_{ij}}=frac{u_{ij}^{p}h_{ij}^{q}}{sumlimits_{k=1}^{c}{u_{kj}^{p}h_{kj}^{q}}} u′ij=k=1∑cukjphkjquijphijq；
- 计算目标函数 J J J，计算公式为 J = ∑ j = 1 N ∑ i = 1 c u i j m ∥ x j − c i ∥ 2 J=sumlimits_{j=1}^{N}{sumlimits_{i=1}^{c}{u_{ij}^{m}{{left| {{x}_{j}}-{{c}_{i}} right|}^{2}}}} J=j=1∑Ni=1∑cuijm∥xj−ci∥2；
- 判断本次迭代与前一次迭代的目标函数的差值，若小于最小误差则退出迭代，若大于误差则返回步骤5；
- 将像素分配给隶属度最大的标签中
%% sFCM 论文代码复现
% Title:Fuzzy c-means clustering with spatial information for image segmentation
% Authors:Keh-Shih Chuang Hong-Long Tzeng Sharon Chen Jay Wu Tzong-Jer Chen
% Computerized Medical Imaging and Graphics
% https://doi.org/10.1016/j.compmedimag.2005.10.001
%% 初始化
clc
clear all
close all
%% 参数设置
error=0.001; % 最小误差
cluster_num=3; % 聚类数
max_iter=100; % 最大迭代次数
m=2; % 隶属度指数
k=5; % 窗口大小
%% 读取图像
f_uint8=imread('testimage.jpg');
[row,col,depth]=size(f_uint8);
N=row*col;
%% 彩色图像转换为灰度图像，并归一化到0到1，便于使用imnoise函数加图像噪声
if depth>1
f=rgb2gray(f_uint8);
f=double(f)/255;
else
f=double(f_uint8)/255;
end
%% 加噪声
f = imnoise(f,'gaussian',0,0.01);
figure, imshow(f),title('噪声图像');
f=f*255;
%% 重组原始图像 便于后续处理
f_reorganize=repmat(f,[1 1 cluster_num]);
%% 初始化隶属度矩阵
U=rand(row,col,cluster_num);
U=U./repmat(sum(U,3),[1 1 cluster_num]);
%% 开始聚类
for iter=1:max_iter
U_m=U.^m;
% 公式（3），更新聚类中心
center(:,iter)=sum(sum(f_reorganize.*U_m),2)./(sum(sum(U_m),2));
% 公式（2）分母部分
d=(f_reorganize-repmat(permute(center(:,iter),[3 2 1]),[row col 1])).^2;
% 公式（2）整体，更新旧的隶属度矩阵
molecular=d.^(1/(m-1));
U=(molecular.*repmat(sum(1./molecular,3),[1 1 cluster_num])).^(-1);
U_m=U.^m;
% 为空间函数 h 分配内存
h=zeros(row,col,cluster_num);
% 扩展隶属度矩阵，便于后续的局部信息操作
U_padded=padarray(U,[k k],'symmetric');
% 公式（3），计算局部隶属度值
for i=-k:k
for j=-k:k
if i==0 && j==0
continue;
end
h=h+U_padded(i+k+1:end+i-k,j+k+1:end+j-k,:);
end
end
% 公式（4），新的隶属度矩阵
U=(U_m.*h)./repmat(sum(U_m.*h,3),[1 1 cluster_num]);
U_m=U.^m;
% 公式（1），更新目标函数
J(iter)=sum(sum(sum(U_m.*d)));
fprintf('第%d次聚类，J=%fn',iter,J(iter));
%判断迭代停止条件
if iter>1 && abs(J(iter)-J(iter-1)) 1 && iter == max_iter && abs(J(iter) - J(iter - 1)) > error
fprintf('目标函数未最小化，提前到达了迭代上限n');
break;
end
end
%% 分割结果显示
[~,cluster_indice]=max(U,[],3);
FCM_result=Label_image(f_uint8,cluster_indice);
figure,imshow(FCM_result),title('分割结果');
J=gather(J);
figure,plot(1:iter,J(1:iter)),title('目标函数曲线');
%%
Vpc = sum(sum(sum(U .^ 2))) / N * 100;
Vpe = -sum(sum(sum(U .* log(U)))) / N * 100;
fprintf('Fuzzy partition coefficient Vpc = %.2f%%n', Vpc);
fprintf('Fuzzy partition entropy Vpe = %.2f%%n', Vpe);
运行结果：
原始图像：
分割结果：（噪声类型为高斯噪声）
目标函数曲线：
V
PC
{{text{V}}_{text{PC}}}
VPC和
V
PE
{{text{V}}_{text{PE}}}
VPE的值：
