GrabCut算法详解：从GMM模型说起

最近正在学习CG，争取有时间就看点论文写写代码。

GrabCut在OpenCv里面是有内置函数的，不过我还是用C++纯手工实现了一边原汁原味的论文，github链接GrabCut，里面有具体的实现代码和步骤，还有实现结果，欢迎来hack

简介

GrabCut是04年的一篇文章，是少数用传统方法~~（对，就是没用神经网络的意思）~~解决前景和背景分离的算法。

需要的前置知识有：

1. 混合多高斯模型

2. 最小割集合划分

下面先介绍这两个模型

备注：本文不含任何数学证明，仅仅介绍算法，如果想要严格证明请移步百度

GMM模型

多维单高斯模型和混合模型

平平无奇的一手公式：维数是DDD，均值是μ\muμ，协方差Σ\SigmaΣ，参数集合θ={μ,Σ}\theta=\{\mu, \Sigma\}θ={μ,Σ}

X∼N(μ,Σ):P(x∣θ)=ϕ(x∣μ,Σ)=1(2π)D/2sqrt∣Σ∣exp{(x−μ)TΣ−1(x−μ)}X\sim N(\mu,\Sigma):P(x|\theta)=\phi(x|\mu,\Sigma)=\frac{1}{(2\pi)^{D/2} sqrt{|\Sigma}|}exp\{(x-\mu)^T\Sigma^{-1}(x-\mu)\} X∼N(μ,Σ):P(x∣θ)=ϕ(x∣μ,Σ)=(2π)D/2sqrt∣Σ∣1​exp{(x−μ)TΣ−1(x−μ)}

其实混合模型就是字面意思，把很多个模型通过一定的概率混合起来，我们首先定义模型的个数KKK，那么参数集合就变成了θ={μk,Σk,πk}\theta=\{\mu_k,\Sigma_k,\pi_k\}θ={μk​,Σk​,πk​}，其中πk\pi_kπk​是第kkk个模型的混合加权系数

P(x∣θ)=∑kKπkϕ(x∣μk,Σk)P(x|\theta)=\sum_k^K\pi_k\phi(x|\mu_k,\Sigma_k) P(x∣θ)=k∑K​πk​ϕ(x∣μk​,Σk​)

混合高斯模型的应用和EM算法

应用很简单：聚类。也就是如果我们假设一堆数据恰好分成K堆并且可以看成每类都是正太分布，那么如果能用GMM模型去拟合这些数据，我们只需要算出每个数据属于某个高斯模型的概率，就可以把这些数据分类了。

拟合的算法有一个经典的EM算法。本质似乎是最大似然概率优化。不过本文通过Local/Global的视角去看待这个算法。

EM算法

输入：数据z1⋯znz_1\cdots z_nz1​⋯zn​和指定的参数KKK

输出：拟合参数θ={μk,Σk,πk}\theta = \{\mu_k,\Sigma_k,\pi_k\}θ={μk​,Σk​,πk​}

1. 初始化参数，包括正则化数据，随机化均值，单位化方差和混合比

2. Esitimate步骤，计算后验概率γjk\gamma_{jk}γjk​

• γjk=πkϕ(zj∣θk)∑kKπkϕ(zi∣θk)\gamma_{jk}=\frac{\pi_k\phi(z_j|\theta_k)} {\sum_k^K\pi_k \phi(z_i|\theta_k)}γjk​=∑kK​πk​ϕ(zi​∣θk​)πk​ϕ(zj​∣θk​)​

3. Maximize步骤，根据后验概率计算新的均值和协方差矩阵以及混合比

• μknew=∑iNγikxi∑iNγik\mu_{k}^{new}=\frac{\sum_i^N\gamma_{ik}x_i}{\sum_i^N \gamma_{ik}}μknew​=∑iN​γik​∑iN​γik​xi​​

• Σknew=∑iNγik(xi−μknew)(xi−μknew)T∑iNγik\Sigma_{k}^{new}=\frac{\sum_i^N\gamma_{ik}(x_i-\mu_k^{new})(x_i-\mu_k^{new} )^T}{\sum_i^N\gamma_{ik}}Σknew​=∑iN​γik​∑iN​γik​(xi​−μknew​)(xi​−μknew​)T​

• πknew=∑iNγik/n\pi_k^{new}=\sum_i^N\gamma_{ik}/nπknew​=∑iN​γik​/n

4. 重复第2步，一直到收敛

这个算法用Local/Global理解非常make sense。也就是先固定参数，计算出对应的每个数据属于某个模型的概率（Local Phase），计算好概率之后，再用对应的概率和数据去更新均值和协方差矩阵以及混合概率。

最小割集合划分

输入：大小为nnn的集合SSS，对于集合的一个划分A/BA/BA/B，定义能量贡献如下：

1. w(i)=[i∈A]Ai+[i∈B]Biw(i)=[i\in A]A_i+[i\in B]B_iw(i)=[i∈A]Ai​+[i∈B]Bi​，也就是每个元素分到某个集合有一定的代价。

2. e(i,j)=[αi≠αj]wije(i,j)=[\alpha_i\neq \alpha_j]w_{ij}e(i,j)=[αi​=αj​]wij​，α\alphaα表示iii的归属。也就是若干对元素，如果他们不在同一个集合，那么有一定的代价

输出：一个划分，最小化能量。

算法：建图，Src→i:Ai,i→Sink:Bi,undirected_edge(i,j,wij)Src\to i:A_i,i\to Sink:B_i,undirected\_edge(i,j,w_{ij})Src→i:Ai​,i→Sink:Bi​,undirected_edge(i,j,wij​)，求最小割即可。

经典算法了，不多解释

GrabCut算法

输入数据：

1. 一张彩色图片，可以转化为三通道向量的二维表数据z1⋯zNz_1\cdots z_Nz1​⋯zN​

2. 用户交互，先圈定一段区域TuT_uTu​，表示这个区域内可能是前景，剩下的部分TB=TˉuT_B=\bar T_uTB​=Tˉu​，表示剩余部分一定是背景

输出数据：

一个αij∈{0,1}\alpha_{ij}\in\{0,1\}αij​∈{0,1}透明度二维表，0表示背景，1表示前景

形式化问题

事实上，上述算法描述是很有问题，因为没有对一个分割做出评估，我并不知道哪个α\alphaα表是最优的。

GrabCut算法的核心就是，把背景和前景分别看成GMM混合模型。并以GMM混合模型参数的拟合程度作为优化目标。

所以优化的参数列表就是，GMM的参数θ={μ(α,k),π(α,k),Σ(α,k)}\theta=\{\mu(\alpha, k),\pi(\alpha, k), \Sigma(\alpha, k)\}θ={μ(α,k),π(α,k),Σ(α,k)}，以及每个像素的透明度αn∈{0,1}\alpha_n\in\{0,1\}αn​∈{0,1}和分类kn∈{0⋯K−1}k_n\in\{0\cdots K-1\}kn​∈{0⋯K−1}（属于第几个高斯模型），也就是每一类的混合高斯模型参数。其中KKK取5。

对于这些参数，我们可以首先定义一个能量函数

U(α,θ,k,z)=∑nD(αn,kn,θ,zn)U(\alpha, \theta, k, z)=\sum_n D(\alpha_n,k_n,\theta,z_n) U(α,θ,k,z)=n∑​D(αn​,kn​,θ,zn​)

其中

D(αn,kn,θ,zn)=−log⁡p(zn∣αn,kn,θ)−log⁡π(αn,kn)=−log⁡π(αn,kn)+log⁡det(Σ(αn,kn))/2+[zn−μ(αn,kn)]TΣ(αn,kn)[zn−μ(αn,kn)]/2D(\alpha_n, k_n,\theta, z_n)=-\log p(z_n|\alpha_n,k_n,\theta)-\log \pi(\alpha_n,k_n)\\ =-\log \pi(\alpha_n,k_n)+\log det(\Sigma(\alpha_n,k_n))/2 +[z_n-\mu(\alpha_n,k_n)]^T\Sigma(\alpha_n,k_n)[z_n-\mu(\alpha_n,k_n)]/2 D(αn​,kn​,θ,zn​)=−logp(zn​∣αn​,kn​,θ)−logπ(αn​,kn​)=−logπ(αn​,kn​)+logdet(Σ(αn​,kn​))/2+[zn​−μ(αn​,kn​)]TΣ(αn​,kn​)[zn​−μ(αn​,kn​)]/2

也就是第k类的混合加权系数以及在这一类的高斯分布的概率密度的log值

当然，除了这个部分以外，分割还关心的是边界的情况，所以论文还定义了另一个能量函数来描述边界的情况

V(α,z)=γ∑(m,n)∈C[αn≠αm]exp−β∣∣zm−zn∣∣2V(\alpha,z)=\gamma \sum_{(m,n)\in C}[\alpha_n\neq \alpha_m]exp-\beta||z_m-z_n||^2 V(α,z)=γ(m,n)∈C∑​[αn​=αm​]exp−β∣∣zm​−zn​∣∣2

CCC是8相邻的像素，γ\gammaγ一般取50，后面那一项实际上表示，边界相差越大，分割得越好，能量越小。

所以β\betaβ参数就应该是图上相邻像素相差的平均情况，也就是β=(2<∣∣zm−zn∣∣2>)−1,<⋅>\beta=(2<||z_m-z_n||^2>)^{-1},<\cdot>β=(2<∣∣zm​−zn​∣∣2>)−1,<⋅>表示在全图上的期望

那么最终的能量优化目标就是

E(αn,kn,θ,zn)=U(α,θ,k,z)+V(α,z)E(\alpha_n,k_n,\theta,z_n)=U(\alpha,\theta, k, z)+V(\alpha, z) E(αn​,kn​,θ,zn​)=U(α,θ,k,z)+V(α,z)

我们先上算法，目前还未实现border matting，只实现了硬分割

算法流程

初始化

对于用户给予的初始分割，我们可以把两个区域内的像素看成数据，并应用GMM算法估计出两个区域内的GMM模型

经典的Local/Global模型

接下来的步骤是不断迭代如下过程，直到能量函数收敛

1. 先固定GMM参数θ\thetaθ，我们算出每个点属于哪个高斯模型的概率最高。

2. 得到了k，我们根据每一类的数据，反推出每一类的参数θ\thetaθ，也就是分别计算每类的均值μ\muμ和方差Σ\SigmaΣ

3. 现在我们确定优化了k,μ,Σk,\mu,\Sigmak,μ,Σ，固定前面这些参数，我们优化α\alphaα以最小化EEE

具体地，我们可以发现，EEE的模型在k,μ,Σk,\mu, \Sigmak,μ,Σ确定的情况下，α=0,1\alpha=0,1α=0,1的代价是固定的，而VVV则相当于是在划分不同类的时候产生贡献。这就是我们的最小割模型，因此得以解决。