GrabCut与BorderMatting的C++实现

最新推荐文章于 2023-02-21 21:17:50 发布

蹦蹦跳跳小米粒

最新推荐文章于 2023-02-21 21:17:50 发布

阅读量3.1k

点赞数 8

CC 4.0 BY-SA版权

分类专栏：计算机视觉 C++ GMM 文章标签： GrabCut MaxFlow OpenCV 计算机视觉

本文链接：https://blog.youkuaiyun.com/weixin_41319239/article/details/91492277

本实验实现了论文《“GrabCut” — Interactive Foreground Extraction using Iterated Graph Cuts》中的图像分割功能，利用迭代的graph cuts进行交互式的前景提取，是对静态图像高效的、交互的前景或背景分割，对图像编辑具有重大的现实意义。

通过在想要提取的物体外围拉出一个包围物体的矩形框，就可以对图像实现较好的分割，提取出我们像要得到的图像部分。论文示例如下：

当前景和背景的颜色比较相近的时候，相似的部分可能会难以分割出来，这时我们可以通过前景或背景的画刷进行交互，通过不同的按键告诉程序自己希望图像该区域所属的部分，并再次进行迭代，以达到更加理想的效果：

除此之外，为了使分割出来得到的图像的分割边缘更加的自然平滑，论文中提出了Border Matting算法，来对得到的图像边缘进行平滑，以得到更加理想的图像效果。论文中的Matting算法对比示例图像如下：

在本项目中，我实现了论文中的GrabCut算法与Border Matting算法，以助教提供的OpenCV框架代码为基础，实现了GrabCut与Border Matting算法的后台代码，通过与OpenCV中相同的方法进行图像的分离操作。

理论分析
1. 图像分割算法概述

图像分割就是将图像分割成若干个特定的、具有独特性质的区域并提出感兴趣目标的技术和过程，是图像处理到图像分析的重要步骤。图像分割较为常用的特征是图像的灰度、颜色、纹理、形状等，这些特征在同一区域呈现出相似性，而在不同区域呈现出明显的差异性。下面首先介绍一下现有的图像分割算法。

2.1.1 基于阈值的图像分割算法

基于阈值的分割算法的基本思想是基于图像的灰度特征来计算一个或者多个灰度阈值，并将图像中每个像素的灰度值与阈值进行比较，最后将比较得到的结果分到合适的类别，因此在本算法中，最重要的部分就是如何找到一个函数来求解最佳的灰度阈值，来保证图像的分割效果。

2.1.2 基于边缘的图像分割算法

基于边缘的分割算法涉及到了图像的边缘提取，是建立在边缘灰度值会呈现出阶跃型或屋顶型变化这一观测基础上的算法。图像的边缘指的是图像中不同区域的边界线的像素点集合，体现出了图像灰度、颜色、纹理等特征的突变，通常情况下基于边缘的图像分割算法指的是基于灰度值的边缘检测，通常体现的特性是灰度值的阶跃性或屋顶性。阶跃性指的是边缘两边的像素的灰度值出现明显的差异，而屋顶性指的是灰度值出现上升或下降的转折处，基于这一特性，可以使用微分算子进行边缘检测来确定边缘，常用的算子有Canny算子、Sobel算子等。

2.1.3 基于区域的图像分割算法

基于区域的分割算法按照图像的相似性准则分为不同的区域，一般有种子区域生长法，区域分裂合并法，分水岭法。

种子区域生长法从一组代表不同图像区域的种子像素开始，不断地将当前像素邻域中符合要求的像素加入到区域像素中，并以该像素为种子继续迭代，直到找不到新像素为止。

区域分裂合并法的基本思想是首先将图像任意分为互不相交的区域，再按照相关准则进行分裂或合并，这种方法既适合灰度图像分割也适用于纹理图像分割。

分水岭法是基于拓扑理论的数学形态学的分割方法，其基本思想是将图像看成地质学上的拓扑打鸡毛，图像上每一点的像素灰度值就是该点的海拔高度，每个局部极小值和其影响的区域为集水盆，集水盆的边缘为分水岭。分水岭算法对边缘微弱的图像有着较好的效果，但是图像中的噪声会使得分水岭算法产生过分割的现象。

2.1.4 基于图论的分割方法

这一类的算法将图像的分割与图像的最小割问题相关链，首先将图像映射为带权的无向图G=<V, E>，图中的每个节点对应于图像中的每一个像素，每条边连接着一对相邻的像素，边的权值表示了像素之间灰度、颜色、纹理等内容的非负相似度，而对图像的一个分割就是对图像的一次裁剪，被分割的每个区域对应着图像的一个子图，分割的最优原则就是让分割后的子图在子图的内部维持相似度最大，而子图之间的相似度最小，本次项目所实现的GrabCut就是基于图论的分割方法。

2.1.5 基于能量泛函的分割方法

这一类的方法主要指的是活动轮廓模型以及在其基础上发展出来的算法，基本思想是使用连续的曲线来表达目标边缘，并且定义一个能量泛函来使得自变量包括边缘曲线，因此分割过程就会变成求解能量泛函最小值的过程，一般可以通过求解函数对应的欧拉方程来实现，当能量达到最小的时候就是目标轮廓的位置。

2.GrabCut算法概述

2.2.1 Graph Cut图割模型

Graph cut图割模型将图像的分割问题与图像的最小割(min cut)问题相结合，将图像转化为一个无向图来表示，无向图为G=<V,E>，V为无向图的顶点，E为无向图的边。其中无向图的顶点和边分为两类，第一类的每一个顶点对应于图像中的每一个像素，而第一类的边为每两个相邻像素之间的边，这种边叫做n-links；而第二类的顶点有两个，叫做S（源点source）点和T（汇点Sink）点，而每个第一类顶点都和这两个终端顶点之间有连接，组成了第二类的边，叫做t-links。

上图即为图像对应的无向图，每个像素在无向图中对应着一个顶点，另外还会有S点和T点，而在图像分割中，我们通常使用s代表图像的前景目标，使用t代表图像的背景目标。图像中的每一条边都有一个非负的权值，每一次图像的切分就意味着消耗掉所切割边的所有权值之和，如果一次切割所删除的边的权值之和最小，则这一次切割就成为最小割，而福特-富克森定理表明，网络的最大流max flow与最小割min cut相等。所以由Boykov和Kolmogorov发明的max-flow/min-cut算法就可以用来获得s-t图的最小割。这个最小割把图的顶点划分为两个不相交的子集S和T，其中s ∈S，t∈ T和S∪T=V 。这两个子集就对应于图像的前景像素集和背景像素集，那就相当于完成了图像分割，所以图像分割的问题就可以转化为无向图中每一条边的权值的确定问题。

图像的分割可以看成是像素的标记问题，目标的标记与背景的标记不同，一般来讲我们将背景标记为0，而将目标标记为1，而我们想要的分割就是要将背景与目标分割开来，而此时分割的cost也是最小的。我们假设每个像素的label为L={l1,l2,...lp}，其中li为0或者1，代表着背景或者前景，当图像的分割为L的时候，图像的能量为：E(L)=aR(L)+B(L)，其中R(L)为区域项，B(L)为边界项，E(L)就是图像的能量函数，我们要求得的图割的最小值来使得能量函数最小。

下面我们分别对区域项和边界项进行分析和解读。

区域项：

其中Rp(Lp)为像素p分配到标签lp的惩罚，在GraphCut中，Rp(Lp)的权值通过比较像素p的灰度值和给定目标的和前景的灰度直方图获得，我们希望像素p分配的标签是其概率最大的标签，但是这时我们希望像素的能量最小，所以我们取概率的负对数。当像素p的灰度值属于前景的概率大于背景的概率时，Rp(1)就会小于Rp(0)，所以像素属于前景的能量会相对更小，因此像素就会被划分到前景中。

边界项：

p和q为邻域像素，B<p,q>为像素p和q之间不连续的惩罚，如果p和q差异非常大，则B<p,q>就接近于0，反之他们属于同一个分类的可能性就比较大，B<p,q>就越大，即能量越大。

综上所述，无向图中的像素顶点之间的边的权值由边界平滑项决定，其权值取决于两个像素之间的差异，另一类边为像素定点与S点和T点的边的权值由区域能量项决定，这时就可以根据mincut对图进行切割，找到使能量最小的图割方案。

2.2.2 GrabCut算法分析

与Graph Cut相比，GrabCut使用了BGR三通道的高斯混合模型进行建模，GrabCut的能量最小化是通过不断进行分割估计和参数模型学习的交互迭代过程，并且允许不完全的标注。

2.2.2.1 RGB颜色空间与GMM建模

在GrabCut中，我们使用5个高斯分量的GMM模型来对前景和背景进行建模，因此存在一个额外的向量k={k1,...kn,...kn}，kn为第n个像素对应于哪个高斯分量，对于每一个像素，可以属于前景高斯分量或者背景高斯分量。整个图像的Gibbs能量为：

其格式与GraphCut的能量方程较为类似，其中U为区域能量项，表示了一个像素被分为前景或者背景像素的惩罚，同样使用像素属于前景或者背景的概率的负对数作为能量值。将GMM模型应用到公式中后，(7)式就变为了(9)式的格式，参数θ有三个，如(10)式所示，分别为每个高斯分量的权重，均值和协方差。所以确定了这三个值后，我们就可以将RGB颜色值带入到前景和背景的GMM求得其属于前景或者背景的概率，即确定了区域能量项。

V为边界能量项，其计算公式为

边界能量项的作用与GraphCut的作用相似，体现了邻域像素m与n之间的不连续惩罚，如果两个像素的差别大，则它们不太可能同属于同一区域，则之间的能量较小，反之则可能同属于同一区域，之间的能量就会较大。在RGB空间中，我们使用欧式距离来衡量两个像素的相似性，即||zm-zn||。其中β由图像的对比度决定，如果图像的对比度较高，则即使两个像素同属一个区域，其欧式距离仍然可能比较大，这时候我们就需要一个β来缩小这种差别，反之亦然，这就使得边界能量项的计算不会在受到图像对比度的影响。γ为50。

2.2.2.2 迭代能量最小化分割图像

GrabCut算法使用了迭代最小化，每一次迭代的过程都使得对前景和背景的GMM参数更优，也就让图像的分割更优。首先用户通过手动框选一个初始化的方框，方框内部的像素为图像的”可能前景”Tu，而方框外部的像素全部被设置为“背景”Tb。我们得到初步的分类后，就可以使用k-means算法将同属于前景和背景的像素聚类为K类，即K个高斯模型，这样每个高斯模型就有了像素样本集，我们就可以根据这些像素求得每一个高斯模型的均值、权值和协方差。

在迭代最小化的过程中，为每个像素分配GMM中的高斯分量，找到每一个像素对应的高斯分量中概率最大的，并将其分配给相应的高斯分量；对于给定的图像，我们需要学习其GMM参数，根据每个高斯分量中已经拥有的像素点，估计其均值和协方差和权重；然后我们对图像进行最大流分割，重复上述操作直到收敛即可。

2.2.2.3 用户交互操作

在经过了一次框选之后，我们可以获得一个初步的分割结果，但是对于图像前景与背景差异不太明显的图像，我们往往难以仅仅通过一次框选就得到令人满意的结果，因此我们需要对图像进行进一步的交互，手动指出图像的部分前景部分和部分背景部分，来帮助程序分析得到更加合适的模型。

3. Border Matting算法概述

我们从前一步的GrabCut中得到的图像前景信息是硬边界的，而Border Matting算法通过“边界消光”机制，使得硬边界分界周围的一个条带上允许透明，以达到让图像边界更加自然的目的。

在Border Matting算法中，我们首先需要找到的是图像的边界信息，因为我们要处理的部分就是图像的边界区域，让图像的过渡更加自然。我们一般选用图像的mask来使用Canny算子进行边缘检测。获取到保存了边缘信息的图像后，我们对图像进行轮廓参数化，找到同属于一个连续轮廓上的像素点，并且保存所有轮廓像素点的信息，为每个像素点分配一个单独的编号，因为后续需要在每一个轮廓上像素点的基础上进行扩展，以找到整个需要处理的条带。在条带查找的阶段，我们一般使用宽度搜索的方式，获取到距离每个中心像素点的一定距离以内的像素点，共同组成了图像前景的边界条带。

为了计算出图像中像素点的alpha值，我们需要确定的参数为sigma与delta，而这两个参数的获取我们使用了动态规划算法，以能量最小化为原则，找到使得能量函数值最小的每个像素点的delta与sigma，最终使用delta与sigma进行计算，就可以得到像素点的alpha值。

能量函数如图所示：

其中分为D项与V项，D项为数据项，V项为平滑项。其中delta与sigma的参数通过最小化该能量函数来获取。

V项的计算方法如下：

其中lambda1取值为50，lambda2取值为1000。

V项的作用是使得alpha的值随着t的增加而平滑的变化，以达到图像的边缘逐渐模糊的效果。在我们使用动态规划算法进行最优的参数值计算的时候，我们通常将delta离散化为30个等级，将sigma离散化为10个等级，在计算中我们使用循环将每一个sigma与delta都代入到轮廓像素中进行计算，直到最后找到使得能量函数取值最小的delta和sigma值并保存。

D项的计算方法如下：

其中N代表了高斯概率密度函数，其中的计算需要使用到均值和协方差，而u和Σ的计算如下：

我们的alpha计算是通过Sigmoid函数进行计算的，其中每个像素使用了其delta和sigma值，获取到整张图像的alphaMask后，将其与原图像进行融合，就得到了border matting后的结果。

3. 代码细节

3.1 GrabCut算法细节

3.1.1 Mask初始化

首先我们创建一个与图像大小相同的mask，且mask中每一个像素点的值为GC_BGD，即均设置为背景。接下来我们判断我们所绘制的矩形的位置，进行放越界限制，并且将矩形内部的像素设置为GC_PR_FGD，即为待定的前景。

3.1.2 初始化GMM模型

我们首先遍历整张图片的所有像素，对于可能是背景的像素全部作为背景的样本像素，而对于所有可能是前景的像素都作为前景的样本像素。

我们定义矩阵bgdSamplesMat，设置其大小为背景样本数*3，3为图像的颜色通道数。我们在这里使用kmeans聚类的方法，将背景的全部样本像素点进行聚类，通常类数为5个，输出为bgdLabels，保存输出样本集中每一个样本对应的类标签，对于前景像素采用相同的处理方法，获得fgdLabels作为标签。我们在得知每一个像素所属的高斯模型之后，下一步需要估计高斯模型中每一个高斯模型的参数，调用AddSample函数，将前景和背景中的每一个像素点加入到其对应的bgdLabels中保存的GMM高斯模型的编号中。所有的像素添加完毕后，我们调用GMMLearning函数，求解每一个高斯模型的权值、均值和方差。权值系数指的是当前样本像素的个数占全部像素个数的比值，并且计算每一个高斯模型的协方差的逆矩阵与行列式，以便于后续步骤的计算。

3.1.3 求解参数

之后我们需要确定一些后续计算中需要应用到的参数。我们在论文中得到gamma为50，构建无向图中使用的lambda为9*gamma，而beta需要我们进行计算得到，beta为Gibbs能量项中的平滑项，用于当图像的对比度较高或较低时对图像相邻像素之间的能量数值进行调整。我们需要分析整幅图像来确定geta的值。首先遍历整个图像的全部像素，计算每个像素与其相邻像素的欧式距离，再计算每两个相邻像素的欧式距离期望，最后计算完成后按照论文中的公式(5)计算得到beta值。

3.1.4 计算无向图边界能量项权重

得到计算所需要的参数后，我们可以开始计算无向图中的每两个顶点之间的边的权重了，这相当于Gibbs能量中的第二项，即平滑项，由于是无向图，我们需要计算八邻域，而对于每一个顶点，我们需要构建的是四个方向的边的权重。首先构造四个矩阵来存储每个点与其左侧、上方、左上方、右上方的权重，然后我们遍历图像中的全部像素，我们首先求出当前像素点与边所连接的像素点的像素值的差，然后调用论文中公式(4)进行计算，并将结果保存在对应的矩阵中。

3.1.5 迭代最小化求解

在求得平滑项后，我们便可以开始迭代求解了。首先我们进行迭代最小化的第一步，为每个像素分配GMM中所属的高斯模型，我们遍历所有像素，如果像素属于前景，则将该像素分配给前景GMM中概率最大的高斯模型；如果像素属于背景，则将该像素分配给背景GMM中概率最大的高斯模型，并将结果保存在与图像大小相同的矩阵compIdxs中。

在初步将所有像素点分配给高斯模型之后，我们调用LearnGMMs来进行迭代最小化求解的第二部，从每个高斯模型的像素样本集中学习每个高斯模型参数。我们循环每个高斯模型，在每一次循环中，我们遍历每一个图像像素点，如果找到图像的所属高斯模型为当前模型，则将该点的信息添加到对应的前景或背景高斯模型中。全部循环完成后，我们再针对每一个高斯模型计算其权重、均值和方差。

在高斯模型参数学习完成后，我们需要建立无向图。在这里我使用了助教在课程网站上提供的Max-flow/min-cut库，调用其中的接口完成对无向图的创建。我们需要通过计算得到的能量项构建图，图的顶点为像素点，边由两部分组成，第一类边为每个像素点顶点与S点和T点的边，其权值由Gibbs能量项的第一项表示，而另一类边为每个顶点与其邻域顶点相连接的边，这一类的边的权值通过Gibbs能量项的第二项平滑项表示。

在计算过程中，我首先遍历图像的全部像素，对于不确定是前景还是背景的顶点，首先计算第一个能量项，其计算公式为论文的公式(3)，对于已经确定是背景或前景的点，我们可以直接赋值：确定为背景的顶点其与S点的权重为0，与T点的权重为lambda；确定为前景的顶点其与S点的权重为lambda，与T点的权重为0，权重计算完毕后我们直接调用函数add_tweights函数，在图中添加该node与source和和sink点的权值。对于像素顶点，我们要获得之前计算得到的平滑项能量值，并且调用add_edge函数将当前顶点与相邻顶点之间的边的权添加到图中。

最后我们使用最小割最大流算法进行分割估计。我们遍历整个mask，如果当前像素点是不确定为前景或者背景的像素点，我们调用what_segment函数进行判断，如果返回值为1则为背景点，如果返回值为0即为前景点。

接下来我们根据自己的需要将上述迭代最小化求解过程进行迭代即可，最终就可以得到一个图像的前景提取结果。

3.2 Border Matting算法细节

3.2.1 图像的边缘检测

首先我们需要对图像的mask进行边缘检测。我们知道mask中的每一个像素值为0、1、2、3中的一个，分别代表着确定背景、确定前景、待定背景、待定前景。我们将mask与1进行与操作，得到新的图像。接下来我们使用Canny算子对图像进行边缘检测，得到图像的边缘信息。

3.2.2 轮廓信息参数化

我们使用深度优先的方法来对轮廓进行参数计算，这一步的作用是找到图像中全部的处于轮廓上的点，并且将其信息保存起来，并且存储当前图像中所有非连续的轮廓的个数，将每个轮廓上的点进行归类，即将位于同一轮廓线上的所有像素点归为同一类，并为每一个轮廓上的像素点分配一个独一无二的编号，为下面分配条带的处理做准备。我们首先遍历轮廓图像的全部像素点，如果当前像素点没有被处理过并且是轮廓上的点，我们就使用深度优先搜索来处理与当前点相邻的周围的点，并将当前的点传入存储了全部轮廓点的vector：contour中，为后续的处理做准备。在这里我们处理每个当前轮廓像素点周围的8个像素点，如果目标点没有超过图像边界、在轮廓线上并且没有被标记过，就将该点作为新的点，并以之为起点再次进行深度优先的搜索，直到所有的轮廓点都被处理完毕。

3.2.3 构建条带信息

在这一步中，我们需要首先初始化Tu，使用一个无序图来进行存储，key值为目标点的坐标信息，value值为目标点。首先，我们需要将所有轮廓上的点添加到strip中，因此我们遍历contour中的全部点，将所有的点加入到无序图strip中，并且添加到向量tmpPointOnStrip中，便于后续的宽搜。遍历完成后，我们以当前已经添加到strip中的点为基础，进行宽度优先的搜索，寻找符合条件的新的点添加到strip中去。我们遍历当前已经添加到条带中的点，并且限制当前点到其条带中心的距离为6，如果超过宽度的限制则不予处理。接下来我们处理当前点相邻的4个点，分别为上、下、左、右四个方向。如果目标点没有超出图像范围，并且没有被处理过，则将其到条带中心的距离绝对值+1，如果目标点是背景点，则距离为负值，设置目标点在条带中所属区域和当前点区域相同，并将目标点加入到条带中。

3.2.4 能量最小化计算

下面需要使用能量最小化方法来求得参数的值，参数即为论文中指出的sigma与delta。首先我们需要使用目标图像的灰度图，由目标图像转化得到。然后我们遍历在轮廓上全部的像素点，计算当前轮廓上的像素点在40*40的区域上的前背景的均值和方差。根据论文，我们将delta离散化为30个等级，将sigma离散化为10个等级，我们在循环中枚举每一个delta与sigma，使用三维数组energyFuncVals记录每个轮廓上的点在每个delta和sigma时的能量值。根据论文，能量函数由两部分组成，分别为D项与V项。

我们首先调用函数计算D项，根据论文的公式(14)，我们需要计算从一个起始点开始，整个轮廓的数据项之和。首先我们计算起始点的数据项值，然后宽度搜索以起始点为中心点的区域，遍历当前点相邻的点，如果目标点属于中心点的区域，就将该点的数据项值添加到sum中，直到循环结束，返回D项结果。由于能量方程是累加的结果，所以如果当前轮廓点是第一个，则直接将得到的D值作为能量值添加到energyFuncVals数组中；否则我们需要枚举index-1时的delta与sigma，如果当前点与上一点属于同一个轮廓，则可以计算V项。V项的计算由论文中的公式(13)决定。

计算完毕后，我们将新的到的D项和V项累加在index-1时的能量值上，并与当前energyFuncVals[index][d0][s0]相比较，如果当前得到的值更小，就用当前值代替之前的值，并记录当前的delta与sigma，直到循环结束。

循环结束后，我们需要找到整个图像总能量的最小值，即能量最小时的delta和sigma值。我们需要遍历所有的deltaLevels和sigmaLevels，找到minEnergy，并记录当前的delta和sigma。记录总能量的delta和sigma后，我们利用该值，取得每个区域的delta和sigma，并保存。

3.2.5 Mask的Alpha值计算

首先我们要遍历轮廓上的所有的点，使用Sigmoid函数(论文中的图6c)计算alpha值，并将其赋值给alpha图像。接下来我们再次宽度优先搜索所有条带上的点，针对当前顶点相邻的4个顶点，如果目标顶点没有超出图像范围并且没有被处理过，则同样使用Sigmoid函数计算该点的alpha值，并添加到alpha图像中。

在获取到包含了Alpha的mask后，我们为了让结果更加明显，对图像进行一次高斯平滑，使结果的显示范围增加。

3.2.6 图像效果展示

在这一阶段，我们将原图与储存了alpha值的mask相融合，主要使用了将每个通道的矩阵与alphaMask相乘，得到的结果merge后就显示出了border matting的结果。

#include "GMM.h"



GMM::GMM(Mat& _model)
{
	if (_model.empty())
	{
		_model.create(1, modelSize*componentsCount, CV_64FC1);
		_model.setTo(Scalar(0));
	}
	model = _model;
	coefs = model.ptr<double>(0);
	mean = coefs + componentsCount;
	cov = coefs + componentsCount * 4;
	for (int i = 0; i < componentsCount; i++)
	{
		covDeterms.push_back(0);
	}
	for (int i = 0; i < componentsCount; i++)
	{
		CalInverseCovAndDeterm(i);
	}
}


GMM::~GMM()
{
}

// 第三次被GMMLearning调用时报错
void GMM::CalInverseCovAndDeterm(int index)
{
	if (coefs[index] > 0)
	{
		double* tmpCovPtr = cov + index * 9;
		double dtrm;
		Mat m = (Mat_<double>(3, 3) << tmpCovPtr[0], tmpCovPtr[1], tmpCovPtr[2],
			tmpCovPtr[3], tmpCovPtr[4], tmpCovPtr[5],
			tmpCovPtr[6], tmpCovPtr[7], tmpCovPtr[8]);
		dtrm = determinant(m);
		covDeterms[index] = dtrm;
		Mat inv;
		invert(m, inv);
		inverseCovMats[index][0][0] = inv.at<double>(0, 0);
		inverseCovMats[index][0][1] = inv.at<double>(0, 1);
		inverseCovMats[index][0][2] = inv.at<double>(0, 2);
		inverseCovMats[index][1][0] = inv.at<double>(1, 0);
		inverseCovMats[index][1][1] = inv.at<double>(1, 1);
		inverseCovMats[index][1][2] = inv.at<double>(1, 2);
		inverseCovMats[index][2][0] = inv.at<double>(2, 0);
		inverseCovMats[index][2][1] = inv.at<double>(2, 1);
		inverseCovMats[index][2][2] = inv.at<double>(2, 2);
	}
}

void GMM::InitGMMLearning()
{
	totalSampleNum = 0;
	for (int i = 0; i < componentsCount; i++)
	{
		sampleNum[i] = 0;
		for (int j = 0; j < 3; j++)
		{
			sums[i][j] = 0;
			for (int k = 0; k < 3; k++)
			{
				prods[i][j][k] = 0;
			}
		}
	}
}

// 为前景后者背景GMM的第i个高斯模型增加样本像素
// 在sum和prod中增加样本值
void GMM::AddSample(int index, Vec3f color)
{
	for (int i = 0; i < 3; i++)
	{
		sums[index][i] += color[i];
		for (int j = 0; j < 3; j++)
		{
			prods[index][i][j] += color[i] * color[j];
		}
	}
	++sampleNum[index];
	++totalSampleNum;
}

// 每个像素有3个通道，因此均值有3个，协方差有9个
void GMM::GMMLearning()
{
	double variance = 0.01;
	for (int i = 0; i < componentsCount; i++)
	{
		int num = sampleNum[i];
		if (num == 0)
			coefs[i] = 0;
		else
		{
			// 权值系数为当前样本像素个数占全部像素个数的比值
			coefs[i] = (double)(num / totalSampleNum);
			double* mVal = mean + 3 * i;
			double* cVal = cov + 9 * i;
			for (int j = 0; j < 3; j++)
			{
				mVal[j] = sums[i][j] / num;
				for (int k = 0; k < 3; k++)
				{
					cVal[j * 3 + k] = prods[i][j][k] / num - mVal[j] * mVal[k];
				}
			}
			Mat tmp = (Mat_<double>(3, 3) << cVal[0], cVal[1], cVal[2],
											cVal[3], cVal[4], cVal[5],
											cVal[6], cVal[7], cVal[8]);
			double d = determinant(tmp);
			// 如果行列式小于等于0，则增加白噪声，防止退化成没有逆矩阵的行列式
			if (d <= 0)
			{
				cVal[0] += variance;
				cVal[4] += variance;
				cVal[8] += variance;
			}
			CalInverseCovAndDeterm(i);
		}
	}
}

int GMM::GetComponent(const Vec3d color)
{
	int idx = 0;
	double max = 0;
	for (int i = 0; i < componentsCount; i++)
	{
		double tmp = GetProbability(i, color);
		if (tmp > max)
		{
			idx = i;
			max = tmp;
		}
	}
	re