【论文阅读】【代码复现】【Matlab】一种用于遥感影像分类的新型自适应模糊局部信息C均值聚类算法(ADFLICM)

[01] Krinidis S, Chatzis V. A robust fuzzy local information C-means clustering algorithm[J]. IEEE transactions on image processing, 2010, 19(5): 1328-1337.
[02] Gong M, Zhou Z, Ma J. Change detection in synthetic aperture radar images based on image fusion and fuzzy clustering[J]. IEEE Transactions on Image Processing, 2011, 21(4): 2141-2151.
[03] Gong M, Liang Y, Shi J, et al. Fuzzy c-means clustering with local information and kernel metric for image segmentation[J]. IEEE transactions on image processing, 2012, 22(2): 573-584.
[04] Zhang H, Wang Q, Shi W, et al. A novel adaptive fuzzy local information C -means clustering algorithm for remotely sensed imagery classification[J]. IEEE Transactions on Geoscience and Remote Sensing, 2017, 55(9): 5057-5068.
[05] Celik T, Lee H K. Comments on “A robust fuzzy local information c-means clustering algorithm”[J]. IEEE transactions on image processing, 2012, 22(3): 1258-1261.

  模糊局部C均值聚类(Fuzzy Local Information C-Means, FLICM) [01] 是由Stelios Krinidis和Vassilios Chatzis于2010年在IEEE Transactions On Image Processing期刊上提出的一种鲁棒的模糊聚类(Fuzzy C-means, FCM)变体算法，该算法是通过构造新的模糊因子提高传统FCM的鲁棒性。由于其算法设计简单、思路新颖，被国内外的许多学者引用和改进，其中引用量已达到1013，它的主要变体算法包括：重新制定的FLICM(Reformulated FLICM, RFLICM) [02]、核加权的FLICM(Kernel Weighted FLICM, KWFLICM) [03]、自适应模糊局部信息C均值(Adaptive Fuzzy Local Information C-Means, ADFLICM)等。题主主要解读和复现ADFLICM算法 [04]。顺带指出，原始FLICM的论文在公式推导过程中存在错误，导致算法的目标函数无法收敛，在论文 [05]中有详细的正确推导步骤，楼主也复现了正确的算法，其聚类效果和原始算法相差不大，因此在此处略去该算法的描述，将重点放在ADFLICM的解读上。

简写：ADFLICM
题目：A Novel Adaptive Fuzzy Local Information C-Means Clustering Algorithm for Remotely Sensed Imagery Classification
作者： Hua Zhang，Qunming Wang，Wenzhong Shi，Ming Hao
题目：一种用于遥感影像分类的新型自适应模糊局部信息 C 均值聚类算法
期刊：IEEE Transactions on Geoscience and Remote Sensing, 2017, 55(9): 5057-5068
doi: 10.1109/TGRS.2017.2702061

  本文提出了一种新型的自适应模糊局部信息C-means（ADFLICM）聚类方法，通过纳入局部空间和灰度信息约束，用于遥感图像分类。ADFLICM方法可以增强传统的模糊C-means算法，产生同质化的分割，同时减少边缘模糊的伪影。ADFLICM的主要贡献是使用了基于像素空间吸引模型的新的模糊局部相似性度量，它可以自适应地确定相邻像素效应的加权因子，而不需要任何实验设置的参数。每个邻域的加权系数对图像内容是完全自适应的，通过使用新的模糊局部相似性度量，自动实现对噪声不敏感和减少边缘模糊假象以保留图像细节之间的平衡。实验中使用了四种不同类型的图像来检验ADFLICM的性能。实验结果表明，ADFLICM比其他四种方法产生更大的准确性，因此为遥感图像的分类提供了一种有效的聚类算法。

  从遥感图像中提取土地覆盖信息是遥感领域的一个常见话题，通常是通过分类来完成的。当训练数据不可用时，无监督聚类被广泛用于遥感图像的分类[1]。许多聚类算法，如K-means[2]，Expectation-最大[3]、ISODATA[4]、K-近邻[5]、马尔科夫随机场（MRF）[6]、模糊c-means（FCM）[7]，以及它们的变种都被用于无监督的分类。其中，FCM[8], [9]是最常用的方法之一。然而，由于有限的空间分辨率、地面物质的复杂性、干扰的多样性或光谱的变化，传统的FCM经常产生含有盐和胡椒噪声的聚类图。
  最近，许多研究人员将局部空间信息纳入常规FCM，以提高聚类性能[10] [21]。最常用的方法之一是修改传统FCM的目标函数以包括空间约束[10] [16], [22] [26]。例如，Pham[23]提出了一种Robust FCM算法，该算法通过包括空间惩罚项来扩展传统的FCM。Ahmed等人[10]，[22]提出了一种FCM_S方法，在FCM目标函数中引入了空间邻接项。FCM_S的一个缺点是耗费时间。为了降低FCM_S的计算复杂性，Chen和Zhang[24]开发了两个变体，FCM_S1和FCM_S2，以简化邻接项的计算。为了进一步加快聚类过程，开发了增强型FCM[25]和快速泛化FCM[11]。然而，这些扩展的FCM算法间接地在原始图像上执行，或者需要一个关键参数来控制对噪声的鲁棒性和保留图像细节的有效性之间的权衡，而这些参数的选择是困难的[10] [13]。为了克服这些问题，Krinidis和Chatzis[13]提出了模糊局部信息C-means（FLICM）。然而，这种方法在识别类边界像素和保留图像细节方面存在一些弱点[14]。为了产生更稳健的结果，Gong等人[14]提出了一个重新表述的FLICM，它引入了一个局部变异系数来代替空间距离作为局部相似度的测量。Li等人[16]提出了一种基于FLICM的带有边缘和局部信息的FCM，以减少边缘退化。在FLICM和上述增强型FLICM算法中，中心像素的识别受其邻近像素的影响很大，而中心像素自身的特征却没有得到充分考虑，未能充分利用局部窗口中封装的局部信息。因此，它们可能会对重要的结构（如区域边界或边缘）和小斑块产生过度光滑的结果。
  为了解决上述问题，本文提出了一种新型的自适应FLICM（ADFLICM）聚类方法，用于遥感图像分类。在ADFLICM中，定义了一个新的模糊局部相似性度量来代替FCM_S中的固定参数α。新的模糊局部相似性度量Sir具有几个特点和优势。 1） S i r {{S}_{ir}} Sir​使用像素空间吸引模型来描述像素之间的关系；2） S i r {{S}_{ir}} Sir​可以通过局部窗口中中心像素和其相邻像素之间的局部空间和灰度关系来自动确定，并且它对局部图像环境具有自适应性，不需要任何人工或经验选择。3）使用 S i r {{S}_{ir}} Sir​，被聚类的像素同时受到其相邻像素和自身特征的影响，这对于去除噪声时保留区域和小斑块的边缘很有帮助；4） S i r {{S}_{ir}} Sir​使所提出的算法相对独立于噪声类型，使其在没有关于噪声的先验知识的情况下成为聚类的有希望的选择。

  假设一幅图像 X = { x 1 , x 2 , ⋯   , x i , ⋯   , x N } X=left{ {{x}_{1}},{{x}_{2}},cdots ,{{x}_{i}},cdots ,{{x}_{N}} right} X={x1​,x2​,⋯,xi​,⋯,xN​}， x i ∈ R n {{x}_{i}}in {{R}^{n}} xi​∈Rn是 n n n维向量空间的数据集， N N N是特征向量的数量（图像中的像素数）， c c c是聚类的数量（ 2 ≤ c ≤ N 2le cle N 2≤c≤N）。

1. 带有空间约束的模糊聚类(FCM_S)
  为了增强传统FCM的稳健性。Ahmed等人[10]引入了一个新的术语，允许一个像素的标签受到其邻居标签的影响。邻近效应作为一个正则器，将解决方案推向零散的同质化标签。FCM_S的目标函数定义如下：
J m = ∑ i = 1 N ∑ k = 1 c u k i m ∥ x i − v k ∥ 2 + α N R ∑ i = 1 N ∑ k = 1 c u k i m ∑ r ∈ N i ∥ x r − v k ∥ 2 {{J}_{m}}=sumlimits_{i=1}^{N}{sumlimits_{k=1}^{c}{u_{ki}^{m}{{left| {{x}_{i}}-{{v}_{k}} right|}^{2}}}}+frac{alpha }{{{N}_{R}}}sumlimits_{i=1}^{N}{sumlimits_{k=1}^{c}{u_{ki}^{m}sumlimits_{rin {{N}_{i}}}{{{left| {{x}_{r}}-{{v}_{k}} right|}^{2}}}}} Jm​=i=1∑N​k=1∑c​ukim​∥xi​−vk​∥2+NR​α​i=1∑N​k=1∑c​ukim​r∈Ni​∑​∥xr​−vk​∥2
其中， x i {{x}_{i}} xi​是第 i i i个像素的灰度值， v k {{v}_{k}} vk​表示第 k k k个簇的原型值， u k i {{u}_{ki}} uki​表示 x i {{x}_{i}} xi​属于第 k k k个簇的模糊成员的程度， m m m是每个模糊成员的权重指数， N R {{N}_{R}} NR​是它的基数， N i {{N}_{i}} Ni​是第 i i i个像素 x i {{x}_{i}} xi​周围寡居像素的集合，像素 x r ( r ∈ N i ) {{x}_{r}}left( rin {{N}_{i}} right) xr​(r∈Ni​)是属于 N i {{N}_{i}} Ni​的邻域像素。参数$alpha $控制来自邻域的影响。FCM_S的细节可以在[10]中找到。
  FCM_S的一个缺点是邻接项的计算很耗时。为了减少计算负担，Chen和Zhang[24]提出了FCM_S1方法。FCM_S1方法，其中FCM_S中的邻接项被简化。修改后的目标函数写成如下：
J m = ∑ i = 1 N ∑ k = 1 c u k i m ∥ x i − v k ∥ 2 + α ∑ i = 1 N ∑ k = 1 c u k i m ∥ x ˉ i − v k ∥ 2 {{J}_{m}}=sumlimits_{i=1}^{N}{sumlimits_{k=1}^{c}{u_{ki}^{m}{{left| {{x}_{i}}-{{v}_{k}} right|}^{2}}}}+alpha sumlimits_{i=1}^{N}{sumlimits_{k=1}^{c}{u_{ki}^{m}{{left| {{{bar{x}}}_{i}}-{{v}_{k}} right|}^{2}}}} Jm​=i=1∑N​k=1∑c​ukim​∥xi​−vk​∥2+αi=1∑N​k=1∑c​ukim​∥xˉi​−vk​∥2
其中 x ˉ i {{bar{x}}_{i}} xˉi​是 x i {{x}_{i}} xi​周围局部窗口内相邻像素的平均值。然而，FCM_S1可能不适合有脉冲噪声的图像[24]。为此，FCM_S1的一个变种，即FCM_S2，在[24]中提出，用中位数过滤的图像来代替平均值过滤的图像，以提高对盐和胡椒噪声等脉冲噪声的鲁棒性。

2.模糊局部信息C-均值聚类算法
   在FCM_S、FCM_S1和FCM_S2的目标函数中，参数 α alpha α平衡了对噪声的鲁棒性和保留图像细节的有效性。FCM_S2的目标函数中，参数 α alpha α平衡了对噪声的鲁棒性和保留图像细节的有效性，它对最终的聚类性能有至关重要的影响。当没有关于噪声的先验知识时，它的选择很困难。在实践中，它通常是根据经验确定的[24]。此外， α alpha α对整个图像的所有邻居窗口都是固定的，局部灰度或空间信息可能被忽略。此外，在FCM_S1和FCM_S2中，使用过滤后的图像可能会导致原始图像的细节丢失。
   为了克服上述局限性，FLICM[13]引入了一个新的模糊因子 G k i {{G}_{ki}} Gki​作为局部相似度测量，以消除噪声并同时保留图像细节：
G k i = ∑ j ∈ N i 1 d i j + 1 ( 1 − u k j ) m ∥ x j − v k ∥ 2 {{G}_{ki}}=sumlimits_{jin {{N}_{i}}}{frac{1}{{{d}_{ij}}+1}{{left( 1-{{u}_{kj}} right)}^{m}}{{left| {{x}_{j}}-{{v}_{k}} right|}^{2}}} Gki​=j∈Ni​∑​dij​+11​(1−ukj​)m∥xj​−vk​∥2
其中 d i j {{d}_{ij}} dij​是像素 i i i和 j j j之间的空间欧几里得距离。将[{{G}_{ki}}]引入传统的FCM，FLICM的目标函数被描述为：
J m = ∑ i = 1 N ∑ k = 1 c [ u k i m ∥ x i − v k ∥ 2 + G k i ] {{J}_{m}}=sumlimits_{i=1}^{N}{sumlimits_{k=1}^{c}{left[ u_{ki}^{m}{{left| {{x}_{i}}-{{v}_{k}} right|}^{2}}+{{G}_{ki}} right]}} Jm​=i=1∑N​k=1∑c​[ukim​∥xi​−vk​∥2+Gki​]
其中 u k i {{u}_{ki}} uki​和 v k {{v}_{k}} vk​的定义为[13]，然而，FLICM在识别类边界像素方面有局限性，重要的结构（如区域边界或边缘）可能被过度平滑[14]。

   在FCM_S、FCM_S1、FCM_S2和FLICM各自优点和局限性的激励下，本文提出了一种新的ADFLICM聚类算法。在这一节中，A节给出了像素间空间吸引力的定义，B节描述了所提出的局部相似度测量方法，C节提供了ADFLICM的一般框架。
A. 像素空间吸引力模型的定义
  吸引模型被证明在描述图像中像素之间的空间相关性方面是有效的[27], [28]。对于两个像素 i i i和 j j j，它们相对于第 k k k个聚类的吸引力与它们的模糊成员 u k i {{u}_{ki}} uki​和 u k j {{u}_{kj}} ukj​成正比，并与两个像素之间的空间距离的平方成反比。因此，两个像素之间的像素空间吸引力 S A i j ( k ) S{{A}_{ij}}left( k right) SAij​(k)可以描述为：
S A i j ( k ) = u k i × u k j D i j 2 S{{A}_{ij}}left( k right)=frac{{{u}_{ki}}times {{u}_{kj}}}{D_{ij}^{2}} SAij​(k)=Dij2​uki​×ukj​​
其中， D i j {{D}_{ij}} Dij​是像素 i i i和 j j j之间的空间距离，本文选择Chebyshev距离为 D i j {{D}_{ij}} Dij​。
B. 新的局部相似性测量
   所提出的相似性措施旨在解决以下问题。1）在对噪声的不敏感性和保存图像细节的有效性之间提供一个适当的权衡；2）权衡应该是自动确定的，而不需要任何人工参数选择；以及 3）该值应同时根据与中心像素的空间距离和灰度级差异而灵活变化。
   在空间吸引模型的基础上，我们引入了一种新的局部相似性测量方法 S i r {{S}_{ir}} Sir​，它同时包含了局部空间和灰度信息。它被定义为：
S i r = S A i r , i ≠ r ; 0 , i = r {{S}_{ir}}=S{{A}_{ir}},ine r;0,i=r Sir​=SAir​,i​=r;0,i=r
其中第 i i i个像素是局部窗口的中心，第 r r r个像素( r ∈ N i rin {{N}_{i}} r∈Ni​)是落入 N i {{N}_{i}} Ni​的邻域像素。本地窗口的邻域结构定义为：
N i = { r ∈ N ∣ 0 < ( a i − a r ) 2 + ( b i − b r ) 2 ≤ Q } {{N}_{i}}=left{ rin Nleft| 0<{{left( {{a}_{i}}-{{a}_{r}} right)}^{2}}+{{left( {{b}_{i}}-{{b}_{r}} right)}^{2}}le Q right. right} Ni​={r∈N∣∣∣​0<(ai​−ar​)2+(bi​−br​)2≤Q}
其中 ( a i , b i ) left( {{a}_{i}},{{b}_{i}} right) (ai​,bi​)和 ( a r , b r ) left( {{a}_{r}},{{b}_{r}} right) (ar​,br​)分别表示像素 i i i和 r r r的坐标， Q Q Q是一个常数，等于 2 L − 1 {{2}^{L-1}} 2L−1。( L L L是邻域的级别)。图1说明了不同级别的邻域结构。请注意，吸引力只存在于中心像素和它在给定的局部窗口中的相邻像素之间，窗口外的其他像素被认为太远而无法对中心像素产生任何吸引力。像素 i i i和 r r r之间的距离 D i r {{D}_{ir}} Dir​可以定义为切比雪夫距离 D i r = max ⁡ ( ∣ a i − a r ∣ , ∣ b i − b r ∣ ) {{D}_{ir}}=max left( left| {{a}_{i}}-{{a}_{r}} right|,left| {{b}_{i}}-{{b}_{r}} right| right) Dir​=max(∣ai​−ar​∣,∣bi​−br​∣)。此外，其他距离指标，如欧氏距离和曼哈顿距离也可以在我们的算法中有利地采用。值得注意的是，局部窗口的形状不仅限于图1中的形状，其他形状如方形也可以在我们的算法中采用。

   显然，局部相似性测量 S i r {{S}_{ir}} Sir​不涉及任何实验调整的参数（除了局部窗口级别 ），以控制图像噪声去除和图像细节保留之间的权衡。邻近效应的加权系数是由空间吸引力自动决定的。 S i r {{S}_{ir}} Sir​措施的引入使得中心像素从其邻居那里得到的影响根据它们的距离 D i r {{D}_{ir}} Dir​而自适应变化。FCM_S2中的参数 α alpha α是全局性的，是一个常数。此外，基于 S i r {{S}_{ir}} Sir​，加权因子不仅受其邻近像素的影响，也受中心像素的影响。这样，当第 r r r个相邻像素的灰度值接近中心像素 i i i的灰度值时，中心像素应该受到这个相邻像素的很大影响，因此， S i r {{S}_{ir}} Sir​应该很大，反之亦然（例如，第 i i i个像素位于边缘区域）。这有助于减少边缘模糊的假象。
C. ADFLICM的总体框架
   在 S i r {{S}_{ir}} Sir​的基础上，提出了ADFLICM用于无监督的遥感图像分类。它在传统FCM的目标函数中加入了局部空间和灰度信息，以增强片状同质分类的平稳性，同时减少边缘模糊效应，其目标函数为ADFLICM的目标函数被描述为：
J m = ∑ i = 1 N ∑ k = 1 c u k i m [ ∥ x i − v k ∥ 2 + 1 N R ∑ r ∈ N i ( 1 − S i r ) ∥ x r − v k ∥ 2 ] {{J}_{m}}=sumlimits_{i=1}^{N}{sumlimits_{k=1}^{c}{u_{ki}^{m}left[ {{left| {{x}_{i}}-{{v}_{k}} right|}^{2}}+frac{1}{{{N}_{R}}}sumlimits_{rin {{N}_{i}}}{left( 1-{{S}_{ir}} right){{left| {{x}_{r}}-{{v}_{k}} right|}^{2}}} right]}} Jm​=i=1∑N​k=1∑c​ukim​[∥xi​−vk​∥2+NR​1​r∈Ni​∑​(1−Sir​)∥xr​−vk​∥2]
其中， N i {{N}_{i}} Ni​的定义见第三节B。 J m {{J}_{m}} Jm​处于局部最小极值的两个必要条件，相对于 u k i {{u}_{ki}} uki​和 v k {{v}_{k}} vk​而言，得到如下：
v k = ∑ i = 1 N u k i m ( x i + 1 N R ∑ r ∈ N i ( 1 − S i r ) ∗ x r ) ( 1 + 1 N R ∑ r ∈ N i ( 1 − S i r ) ) ∑ i = 1 N u k i m {{v}_{k}}=frac{sumnolimits_{i=1}^{N}{u_{ki}^{m}left( {{x}_{i}}+frac{1}{{{N}_{R}}}sumnolimits_{rin {{N}_{i}}}{left( 1-{{S}_{ir}} right)*{{x}_{r}}} right)}}{left( 1+frac{1}{{{N}_{R}}}sumnolimits_{rin {{N}_{i}}}{left( 1-{{S}_{ir}} right)} right)sumnolimits_{i=1}^{N}{u_{ki}^{m}}} vk​=(1+NR​1​∑r∈Ni​​(1−Sir​))∑i=1N​ukim​∑i=1N​ukim​(xi​+NR​1​∑r∈Ni​​(1−Sir​)∗xr​)​ u k i = 1 ∑ j = 1 c ( ∥ x i − v k ∥ 2 + 1 N R ∑ r ∈ N i ( 1 − S i r ) ∥ x r − v k ∥ 2 ∥ x i − v j ∥ 2 + 1 N R ∑ r ∈ N i ( 1 − S i r ) ∥ x r − v j ∥ 2 ) 1 / ( m − 1 )    {{u}_{ki}}=frac{1}{sumnolimits_{j=1}^{c}{{{left( frac{{{left| {{x}_{i}}-{{v}_{k}} right|}^{2}}+frac{1}{{{N}_{R}}}sumlimits_{rin {{N}_{i}}}{left( 1-{{S}_{ir}} right){{left| {{x}_{r}}-{{v}_{k}} right|}^{2}}}}{{{left| {{x}_{i}}-{{v}_{j}} right|}^{2}}+frac{1}{{{N}_{R}}}sumlimits_{rin {{N}_{i}}}{left( 1-{{S}_{ir}} right){{left| {{x}_{r}}-{{v}_{j}} right|}^{2}}}} right)}^{{1}/{left( m-1 right)};}}}} uki​=∑j=1c​(∥xi​−vj​∥2+NR​1​r∈Ni​∑​(1−Sir​)∥xr​−vj​∥2∥xi​−vk​∥2+NR​1​r∈Ni​∑​(1−Sir​)∥xr​−vk​∥2​)1/(m−1)1​

   论文中采用马尔科夫场随机生成的合成图像和真实遥感图像作为实验数据，进行算法的验证。题主这里用三种类型的图像：合成图像、彩色图像、遥感图像和MR图像，用复现的算法测试其鲁棒性（注意，原始代码的测试结果并不理想[如E所示]，因此题主将区域级信息引入目标函数，提高鲁棒性[如A、B、C、D所示]）：
A. 合成图像实验

B.彩色图像实验

B.遥感图像实验

D.MR图像实验

E.原始代码分割结果

   题主根据论文公式以及伪代码，自主复现的原创Matlab代码如下：

%% ADFLICM
% 一种用于遥感影像分类的新型自适应模糊局部信息 C 均值聚类算法
% 作者: Hua Zhang，Qunming Wang，Wenzhong Shi，Ming Hao
% A Novel Adaptive Fuzzy Local Information C-Means Clustering Algorithm
% for Remotely Sensed Imagery Classification
% IEEE Transactions on Geoscience and Remote Sensing, 2017, 55(9): 5057-5068
% doi: 10.1109/TGRS.2017.2702061
%% 初始化
clc
clear all
close all
%% 参数设定
error=0.001;    % 最小误差
cluster_num=4;  % 聚类数
max_iter=100;   % 最大迭代次数
m=2;    % 隶属度指数
k=3;    % 局部窗口
%% 读入图像
f_uint8=imread('2.bmp');
[row,col,depth]=size(f_uint8);
%% 彩色图像转换为灰度图像，并归一化到0到1，便于使用imnoise函数加图像噪声
if depth>1
    f=rgb2gray(f_uint8);
    f=double(f)/255;
else
    f=double(f_uint8)/255;
end
%% 加噪声
density=0.10;
f = imnoise(f,'gaussian',0,density);
f = imnoise(f,'salt & pepper',density);
f = imnoise(f,'speckle',density);
figure, imshow(f),title('噪声图像');
f=f*255;
%% 初始化隶属度矩阵
U=rand(row,col,cluster_num);
U=U./repmat(sum(U,3),[1 1 cluster_num]);
%% 扩展原始图像，便于后续的局部信息操作
f_padded=padarray(f,[k k],'symmetric');
f_padded=repmat(f_padded,[1 1 cluster_num]);
f_dense=repmat(f,[1 1 cluster_num]);
%% 开始聚类
for iter=1:max_iter
    U_m=U.^m; % U的m次方
    % 扩展隶属度矩阵，便于以后局部运算
    U_padded=padarray(U,[k k],'symmetric');
    % 为公式（9）的后半部分开辟空间
    SA1=zeros(row,col,cluster_num);
    SA=zeros(row,col,cluster_num); 
    for i=-k:k
        for j=-k:k
            if i==0 && j==0
                continue;
            end
            % 根据公式（6）计算空间吸引模型S_ir
            S=(1/(sqrt(i*i+j*j)).^2).*(U_padded(i+k+1:end+i-k,j+k+1:end+j-k,:).*U);
            % 计算公式（9）的后半部分
            SA1=SA1+(1-S);
            SA=SA+((1-S).*f_padded(i+k+1:end+i-k,j+k+1:end+j-k,:));
        end
    end
    SA1=SA1/(2*k+1);
    SA=SA/(2*k+1);
    % 计算公式（9），更新聚类中心
    center(:,iter)=sum(sum((f_dense+SA).*U_m),2)./(sum(sum(U_m+U_m.*SA1),2)); 
    % 为公式（8）和（10）的后半部分开辟空间
    G=zeros(row,col,cluster_num);
    for i=-k:k
        for j=-k:k
            if i==0 && j==0
                continue;
            end
            % 计算公式（8）和（10）中的||x_r-v_k||^2，更新局部距离
            d_local=(f_padded(i+k+1:end+i-k,j+k+1:end+j-k,:)-repmat(permute(center(:,iter),[3 2 1]),[row col 1])).^2;
            % 根据公式（6）计算空间吸引模型S_ir
            S1=(1/(sqrt(i*i+j*j)).^2).*(U_padded(i+k+1:end+i-k,j+k+1:end+j-k,:).*U);
            % 计算公式(17)中求和符号里面的内容，更新模糊因子G
            G=G+(1-S1).*d_local;
        end
    end
    G=G/(2*k+1);
    % 计算||x_i-v_k||^2，更新欧式距离
    d=(f_dense-repmat(permute(center(:,iter),[3 2 1]),[row col 1])).^2;
    % 计算公式(10)，更新聚类隶属度
    U_old=U; % 保存旧隶属度矩阵，用于迭代停止条件的计算
    fenzi=(d+G).^(1/(m-1));
    U=(fenzi.*repmat(sum(1./fenzi,3),[1 1 cluster_num])).^(-1);
    U_m=U.^m;
    % 计算公式(18)，更新目标函数
    J(iter)=sum(sum(sum(U_m.*(d+G))));
    fprintf('第%d次聚类，J=%fn',iter,J(iter));
    %判断迭代停止条件
    if iter>1 && abs(max(max(max(U-U_old))))1 && iter==max_iter && abs(max(max(max(U-U_old))))>error
        fprintf('目标函数未最小化，提前到达了迭代上限n');
        break;
    end
end
%% 分割结果显示
[~,cluster_indice]=max(U,[],3);
FCM_result=Label_image(f_uint8,cluster_indice);
figure,imshow(FCM_result);
J=gather(J);
figure,plot(1:iter,J(1:iter));