【图像分割】【FCM】【经典文献解读】基于隐Markov随机场模型的增强型空间约束FCM图像分割 HMRF-FCM

参考文献
Chatzis S. P., Varvarigou T. A. A Fuzzy Clustering Approach Toward Hidden Markov Random Field Models for Enhanced Spatially Constrained Image Segmentation. IEEE Transactions on Fuzzy Systems 16(5) (2008) 1351-1361
该算法首次将隐Markov随机场模型引入FCM，达到了较好的分割效果。

1、隐Markov随机场模型

随机场是什么？打一个很简单的比喻:现在有一个农民，家里有几块地，开春了，要开始种庄稼，每一块地里种什么庄稼，都具有一个概率，所有位置的概率构成了随机场。
什么是Markov随机场？简单地说，就是某一块地种什么庄稼，只与这块地周围的地种什么庄稼有关，与更远的土地没有关系。

下面给出数学表达。对于图像分割，实际上就是给每个像素打标签，标签表明了每个像素被分进了哪个类别。假设想要将10个像素分成3个类，每个像素只能分进1类中。
令$Q$为标签集
Q = { 1 , 2 , 3 } Q=\{1,2,3\} Q={1,2,3}
其中1表示第一个类，2表示第二个类，依次类推。
令S为索引集
S = { 1 , 2 , . . . , 10 } S=\{1,2,...,10\} S={1,2,...,10}
其中1表示第一个像素，2表示第二个像素，依次类推。
则对于任意$j\in S$，有$x_j\in Q$，$x_j$为第$j$个像素所属的类别，可以为1，也可以为2，等等。
此时定义一个概率分布
$$p(x_j)>0$$
表示第$j$个像素为某个类别的概率。该概率分布，就是随机场。
OK，现在讨论Markov随机场~
前面说了，在Markov随机场中，某一个像素属于什么类别只与它周围像素属于什么类别有关系，那么显然，必须首先定义一个能够描述“周围”的方式，在数学上，这种方式为邻域。我们知道，每一个像素都拥有与其紧邻的像素（不考虑位于边缘的像素），即每一个像素都具有邻域。

令$\partial$表示所有像素的邻域，即邻域系统。那么现在，就可以用数学语言来描述Markov随机场
$$p(x_j)=p(x_j|X_{{\partial }_j})$$
其中，$X_{{\partial }_j}$表示一个集合，包含第$j$个像素邻域中所有像素的类别。
然而，图像是二维的，从而，$j$实际上是一个二维坐标，所以接下来，为了更加详细地描述，需要使用联合概率分布。Hammersley–Clifford理论证明，Markov随机场与Gibbs随机场的联合概率分布是等价的，则Markov随机场的联合概率分布为
$$p(x_j)=p(x_j|X_{{\partial }_j})=\frac{\exp (-U(x_j))}{Z}$$
其中，$Z$是一个使$p(x_j)$归一化的常数，$U(x_j)$叫做能量函数。
$$U(x_j)=\sum _C{V}_C(x_j)$$
其中，$V_C(x_j)$表示第$j$个像素邻域中某一个基团C的势能，OK，如何表示势能呢，只能说八仙过海，各显神通，目前提出了多种势能表示方法，OK，那么什么是基团，简单地说，就是邻域里面，一个、两个或三个相邻的像素可以拉帮结伙，构成一个个小集团，分别为一阶基团、二阶基团或三阶基团。
人们发现，在实际应用中，$Z$的计算量较大，需要一种近似计算方法减少计算量。1975年，J.Besag提出一种伪似然估计法，即
$$p(x_j)=p(x_j|X_{{\partial }_j})=\frac{\exp (-{\sum }_{j\in C}{V}_C(x_j))}{{\sum }_{x_j=1}^3\exp (-{\sum }_{j\in C}{V}_C(x_j))}$$
现在，Markov随机场已经被表示了出来。

然而，实际生活经验告诉我们，将所有像素一视同仁随便分类显然不科学，不同的像素需要根据其灰度值来进行分类，OK，现在引入另一个条件，像素灰度值。
令集合$R$表示所有灰度值，有
R = { 0 , 1 , 2 , 3 , 4 , 5 , . . . , 255 } R=\{0,1,2,3,4,5,...,255\} R={0,1,2,3,4,5,...,255}
对于任意$j\in S$，有$y_j\in R$，$y_j$表示第$j$个像素的灰度。对于一张图像，每个像素的灰度，一定是事先就确定的。那么将灰度$y_j$与像素类别$x_j$放在一起，构成的联合概率分布$p(y_j,x_j)$，就是隐Markov随机场。
接下来如何求解$p(y_j,x_j)$，是得到隐Markov随机场的关键。
众所周知
$$p(y_j,x_j)=p(y_j|x_j)p(x_j)$$
其中，$p(x_j)$已经可以表示，它就是Gibbs随机场的联合概率分布，叫做隐含场，$p(y_j)$叫做观测场。需要求解$p(y_j|x_j)$，即已知第$j$个像素的类别，它的灰度是多少的概率。一般地，可以利用正态分布，描述$p(y_j|x_j)$，因为正态分布有参数，所以引入$\theta$表示其所有参数，即
$$p(y_j|x_j)=p(y_j|x_j;{\theta }_{x_j})=N_{x_j}({\mu }_{x_j},{\sigma }_{x_j})$$
其中，${\mu }_{x_j}和{\sigma }_{x_j}$分别表示第$x_j$类中所有像素灰度值的均值与方差，一定记住，$x_j$表示的是类别标签。将$p(y_j|x_j)$画成图像的话，则横坐标为灰度值，纵坐标为概率值，由于分3类，则有3条高斯概率密度曲线。
现在，隐Markov随机场也被表示完成。

OK，现在将隐Markov随机场引入FCM中

2、基于隐Markov随机场的FCM—HMRF-FCM

首先，该算法利用了Ichihashi等人提出的一种基于KL信息惩罚的FCM方法，目标函数如下
$$J=\sum _{i=1}^3\sum _{j=1}^{10}{u}_{ij}{d}_{ij}+\lambda \sum _{i=1}^3\sum _{j=1}^{10}{u}_{ij}\ln (\frac{u_{ij}}{{\pi }_{ij}})$$
其中，
$${\pi }_{ij}=p(x_j)=p(x_j=i|X_{{\partial }_j})=\frac{\exp (-{\sum }_{j\in C}{V}_C(a_{ij}))}{{\sum }_{h=j}^3\exp (-{\sum }_{j\in C}{V}_C(a_{hj}))}$$
$a_{ij}=(x_j=i,X_{{\theta }_j})$，由于$X_{{\theta }_j}$是一个集合，在实际应用中，不容易直接使用，一般利用平均场估计方法得到$X_{{\theta }_j}$的估计，记为$x_{{\theta }_j}^{\prime }$，则$a_{ij}=(x_j=i,x_{{\theta }_j}^{\prime })$
本算法为简化操作，直接利用8邻域中的一阶基团像素标签表达能量函数$U(x_j)$，即
$$U(x_j)=-\beta \sum _{j=1}^{10}\sum _{l=1}^8\delta (x_j-x_l)$$
其中，$x_l$为$x_j$的邻域${\partial }_j$中的像素标签，则有
$${\pi }_{ij}=p(x_j)=p(x_j=i|X_{{\partial }_j})=\frac{\beta \exp (-{\sum }_{l\in {\partial }_j}\delta (i-x_l))}{{\sum }_{h=1}^3\beta \exp (-{\sum }_{l\in {\partial }_j}\delta (h-x_l))}$$
其中$\delta (x_j-x_l)$为克罗内克函数，
$$\delta (x_j,x_l)=\{\begin{array}{l}1,x_j=x_l\\ 0,otherwise\end{array}$$
通过拉格朗日乘子法，对$J$最小化，得到隶属度$r_{ij}$、聚类中心${\mu }_i$、方差${\sigma }_i$的更新公式
$$r_{ij}=\frac{{\pi }_{ij}\exp (-\frac{1}{\lambda }{d}_{ij})}{{\sum }_{h=1}^3{\pi }_{ij}\exp (-\frac{1}{\lambda }{d}_{ij})}，d_{ij}=\frac{w}{2}\ln (2\pi )+\frac{1}{2}\ln (|{\mu }_i|)+\frac{(y_j-{\mu }_i{)}^2}{2{\mu }_i}$$
$${\mu }_i=\frac{{\sum }_{j=1}^{10}{r}_{ij}{y}_j}{{\sum }_{j=1}^{10}{y}_j}$$
$${\sigma }_i=\frac{{\sum }_{j=1}^{10}{r}_{ij}(y_j-{\mu }_i{)}^2}{{\sum }_{j=1}^{10}{r}_{ij}}$$
算法步骤为：
1、随机初始化隶属度$r_{ij}$
2、根据隶属度，得到所有像素的类别标签 x j = arg max ⁡ i { r i j } x_j=\argmax_i\{r_{ij}\} xj​=iargmax​{rij​}
3、计算先验概率${\pi }_{ij}$
4、计算聚类中心${\mu }_i$和方差${\sigma }_i$
4、计算隶属度$r_{ij}$
5、若收敛，则输出$r_{ij}$，否则返回2
最终的分类标签为 x j = arg max ⁡ i { r i j } x_j=\argmax_i\{r_{ij}\} xj​=iargmax​{rij​}