本实验实现了论文《“GrabCut” — Interactive Foreground Extraction using Iterated Graph Cuts》中的图像分割功能，利用迭代的graph cuts进行交互式的前景提取，是对静态图像高效的、交互的前景或背景分割，对图像编辑具有重大的现实意义。
通过在想要提取的物体外围拉出一个包围物体的矩形框，就可以对图像实现较好的分割，提取出我们像要得到的图像部分。论文示例如下：
当前景和背景的颜色比较相近的时候，相似的部分可能会难以分割出来，这时我们可以通过前景或背景的画刷进行交互，通过不同的按键告诉程序自己希望图像该区域所属的部分，并再次进行迭代，以达到更加理想的效果：
除此之外，为了使分割出来得到的图像的分割边缘更加的自然平滑，论文中提出了Border Matting算法，来对得到的图像边缘进行平滑，以得到更加理想的图像效果。论文中的Matting算法对比示例图像如下：
在本项目中，我实现了论文中的GrabCut算法与Border Matting算法，以助教提供的OpenCV框架代码为基础，实现了GrabCut与Border Matting算法的后台代码，通过与OpenCV中相同的方法进行图像的分离操作。
- 理论分析
- 图像分割算法概述
图像分割就是将图像分割成若干个特定的、具有独特性质的区域并提出感兴趣目标的技术和过程，是图像处理到图像分析的重要步骤。图像分割较为常用的特征是图像的灰度、颜色、纹理、形状等，这些特征在同一区域呈现出相似性，而在不同区域呈现出明显的差异性。下面首先介绍一下现有的图像分割算法。
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;
}
}
return idx;
}
double GMM::GetProbability(const Vec3d color)
{
double result = 0;
for (int i = 0; i < componentsCount; i++)
{
result += coefs[i] * GetProbability(i, color);
}
return result;
}
// 计算当前像素与当前高斯分量的均值的差
// 使用diff转置*inverseCovMats*diff求得mul
double GMM::GetProbability(int index, const Vec3d color)
{
double result = 0;
if (coefs[index] > 0)
{
if (covDeterms[index] > 0)
{
double* mVal = mean + 3 * index;
Vec3d diff = color;
for (int i = 0; i < 3; i++)
{
diff[i] -= mVal[i];
}
double mul = 0;
for (int i = 0; i < 3; i++)
{
for (int j = 0; j < 3; j++)
{
mul += diff[i] * diff[j] * inverseCovMats[index][j][i];
}
}
// 多维变量的联合高斯密度函数求解概率
result = 1.0 / sqrt(covDeterms[index]) * exp(-0.5f*mul);
}
else
{
}
}
return result;
}
#pragma once
#include<iostream>
#include<vector>
#include<opencv2\opencv.hpp>
using namespace std;
using namespace cv;
class GMM
{
public:
GMM(Mat& _model);
~GMM();
static const int componentsCount = 5;
void CalInverseCovAndDeterm(int index);
void InitGMMLearning();
void AddSample(int index, Vec3f color);
void GMMLearning();
int GetComponent(const Vec3d color);
double GetProbability(const Vec3d color);
double GetProbability(int index, const Vec3d color);
private:
int modelSize = 13;
Mat model;
double* cov;
double* coefs;
double* mean;
int sampleNum[componentsCount];
int totalSampleNum;
// 储存了每个component的协方差的逆矩阵
double inverseCovMats[componentsCount][3][3] ;
// 储存了每个component的协方差的行列式
vector<double> covDeterms;
// 每个component的sum
double sums[componentsCount][3];
double prods[componentsCount][3][3];
};
#include "GrabCut.h"
void GrabCut2D::CheckRectAndMask(Mat &_mask, Rect rect, int mode, Size imgSize)
{
if (mode == GC_INIT_WITH_RECT)
{
InitMaskWithRect(_mask, rect, imgSize);
}
else if (mode == GC_INIT_WITH_MASK)
{
}
else
{
}
}
void GrabCut2D::InitMaskWithRect(Mat & _mask, Rect rect, Size imgSize)
{
if (_mask.empty())
{
_mask.create(imgSize, CV_8UC1);
_mask.setTo(Scalar(GC_BGD));
}
else
{
_mask.setTo(GC_BGD);
}
// 防止出现矩形越界
if (rect.x < 0)
rect.x = 0;
if (rect.y < 0)
rect.y = 0;
if (rect.width > imgSize.width - rect.x)
rect.width = imgSize.width - rect.x;
if (rect.height > imgSize.height - rect.y)
rect.height = imgSize.height - rect.y;
_mask(rect).setTo(Scalar(GC_PR_FGD));
//cout << GC_PR_FGD << endl;
//cout << Scalar(GC_PR_FGD) << endl;
}
GrabCut2D::~GrabCut2D(void)
{
}
// k-means聚类初始化两个GMM
void GrabCut2D::InitGMMs(Mat& img, Mat& mask, GMM& bgdGMM, GMM& fgdGMM)
{
const int iterCnt = 10;
const int kMeansType = KMEANS_PP_CENTERS;
// 记录前景和背景的像素样本集中每个像素对应GMM中的哪一个高斯模型
Mat bgdLabels, fgdLabels;
vector<Vec3f> bgdSamples, fgdSamples;
// 遍历整个图片的全部像素
for (int i = 0; i < img.rows; i++)
{
for (int j = 0; j < img.cols; j++)
{
// mask中所有可能是背景的像素都作为背景的样本像素
if (mask.at<uchar>(i, j) == GC_BGD || mask.at<uchar>(i, j) == GC_PR_BGD)
{
bgdSamples.push_back((Vec3d)img.at<Vec3b>(i, j));
}
// mask中所有可能是前景的像素都作为前景的样本像素
else
{
fgdSamples.push_back((Vec3d)img.at<Vec3b>(i, j));
}
}
}
Mat bgdSamplesMat;
bgdSamplesMat.create((int)bgdSamples.size(), 3, CV_32FC1);
bgdSamplesMat.setTo(bgdSamples[0][0]);
kmeans(bgdSamplesMat, GMM::componentsCount, bgdLabels,
TermCriteria(CV_TERMCRIT_ITER, iterCnt, 0), 0, kMeansType);
Mat fgdSamplesOutput;
fgdSamplesOutput.create((int)fgdSamples.size(), 3, CV_32FC1);
fgdSamplesOutput.setTo(fgdSamples[0][0]);
kmeans(fgdSamplesOutput, GMM::componentsCount, fgdLabels,
TermCriteria(CV_TERMCRIT_ITER, iterCnt, 0), 0, kMeansType);
// 已经确定每个像素所属的高斯模型，进而估计GMM中每个高斯模型参数
bgdGMM.InitGMMLearning();
fgdGMM.InitGMMLearning();
for (int i = 0; i < bgdSamples.size(); i++)
{
bgdGMM.AddSample(bgdLabels.at<int>(i, 0), bgdSamples[i]);
}
bgdGMM.GMMLearning();
for (int i = 0; i < fgdSamples.size(); i++)
{
fgdGMM.AddSample(fgdLabels.at<int>(i, 0), fgdSamples[i]);
}
fgdGMM.GMMLearning();
//Mat bgdSamplesOutput((int)bgdSamples.size(), 3, CV_32FC1, &bgdSamples[0][0]);
}
double GrabCut2D::CalBeta(const Mat& img)
{
double beta = 0;
for (int i = 0; i < img.rows; i++)
{
for (int j = 0; j < img.cols; j++)
{
// 左侧
if (j > 0)
{
Vec3d diff = (Vec3d)img.at<Vec3b>(i, j) - (Vec3d)img.at<Vec3b>(i, j - 1);
beta += diff.dot(diff);
}
// 左上方
if (i > 0 && j > 0)
{
Vec3d diff = (Vec3d)img.at<Vec3b>(i, j) - (Vec3d)img.at<Vec3b>(i - 1, j - 1);
beta += diff.dot(diff);
}
// 上方
if (i > 0)
{
Vec3d diff = (Vec3d)img.at<Vec3b>(i, j) - (Vec3d)img.at<Vec3b>(i - 1, j);
beta += diff.dot(diff);
}
// 右上方
if (i > 0 && j < img.cols - 1)
{
Vec3d diff = (Vec3d)img.at<Vec3b>(i, j) - (Vec3d)img.at<Vec3b>(i - 1, j + 1);
beta += diff.dot(diff);
}
}
}
if (beta < 0)
{
beta = 0;
}
else
{
beta = 1.0f / (2 * beta / (4 * img.cols * img.rows - 3 * img.cols - 3 * img.rows + 2));
}
return beta;
}
// 计算图中每个像素顶点的权值，即与其相连的8个像素的权值
void GrabCut2D::CalAllWeights(const Mat& img, Mat& leftWeight, Mat& upLeftWeight, Mat& upWeight, Mat& upRightWeight,
double beta, double gamma)
{
const double gammaNear = gamma;
const double gammaFar = gamma / sqrt(2.0);
leftWeight.create(img.rows, img.cols, CV_64FC1);
upLeftWeight.create(img.rows, img.cols, CV_64FC1);
upWeight.create(img.rows, img.cols, CV_64FC1);
upRightWeight.create(img.rows, img.cols, CV_64FC1);
for (int i = 0; i < img.rows; i++)
{
for (int j = 0; j < img.cols; j++)
{
// 左侧权重
if (j > 0)
{
Vec3d diff = (Vec3d)img.at<Vec3b>(i, j) - (Vec3d)img.at<Vec3b>(i, j - 1);
leftWeight.at<double>(i, j) = gammaNear * exp(diff.dot(diff) * (beta * -1));
}
else
leftWeight.at<double>(i, j) = 0;
// 上方权重
if (i > 0)
{
Vec3d diff = (Vec3d)img.at<Vec3b>(i, j) - (Vec3d)img.at<Vec3b>(i - 1, j);
upWeight.at<double>(i, j) = gammaNear * exp(diff.dot(diff) * (beta * -1));
}
else
upWeight.at<double>(i, j) = 0;
// 左上方权重
if (i > 0 && j > 0)
{
Vec3d diff = (Vec3d)img.at<Vec3b>(i, j) - (Vec3d)img.at<Vec3b>(i - 1, j - 1);
upLeftWeight.at<double>(i, j) = gammaFar * exp(diff.dot(diff) * (beta * -1));
}
else
upLeftWeight.at<double>(i, j) = 0;
// 右上方权重
if (j < img.cols - 1 && i > 0)
{
Vec3d diff = (Vec3d)img.at<Vec3b>(i, j) - (Vec3d)img.at<Vec3b>(i - 1, j + 1);
upRightWeight.at<double>(i, j) = gammaFar * exp(diff.dot(diff) * (beta * -1));
}
else
upRightWeight.at<double>(i, j) = 0;
}
}
}
// 迭代最小化的步骤一，给每个像素分配GMM中所属的高斯模型
void GrabCut2D::AssignGMMsComponents(const Mat& img, const Mat& mask, GMM& bgdGMM, GMM& fgdGMM,
Mat& compIdxs)
{
for (int i = 0; i < img.rows; i++)
{
for (int j = 0; j < img.cols; j++)
{
if (mask.at<uchar>(i, j) == GC_BGD || mask.at<uchar>(i, j) == GC_PR_BGD)
{
Vec3d v = (Vec3d)(img.at<Vec3b>(i, j));
compIdxs.at<int>(i, j) = bgdGMM.GetComponent(v);
}
else
{
Vec3d v = (Vec3d)(img.at<Vec3b>(i, j));
compIdxs.at<int>(i, j) = fgdGMM.GetComponent(v);
}
}
}
}
// 迭代最小化算法step 2：从每个高斯模型的像素样本集中学习每个高斯模型的参数
void GrabCut2D::LearnGMMs(const Mat & img, const Mat & mask, const Mat & compIdxs, GMM & bgdGMM, GMM & fgdGMM)
{
bgdGMM.InitGMMLearning();
fgdGMM.InitGMMLearning();
for (int index = 0; index < GMM::componentsCount; index++)
{
for (int i = 0; i < img.rows; i++)
{
for (int j = 0; j < img.cols; j++)
{
if (compIdxs.at<int>(i, j) == index)
{
if (mask.at<uchar>(i, j) == GC_BGD || mask.at<uchar>(i, j) == GC_PR_BGD)
{
bgdGMM.AddSample(index, img.at<Vec3b>(i, j));
}
else
{
fgdGMM.AddSample(index, img.at<Vec3b>(i, j));
}
}
}
}
}
bgdGMM.GMMLearning();
fgdGMM.GMMLearning();
}
void GrabCut2D::ConstructGraph(const Mat & img, const Mat & mask, GMM & bgdGMM, GMM & fgdGMM, double lambda, const Mat & leftWeight, const Mat & upLeftWeight, const Mat & upWeight, const Mat upRightWeight, Graph<double, double, double>& graph)
{
for (int i = 0; i < img.rows; i++)
{
for (int j = 0; j < img.cols; j++)
{
int nodeIndex = graph.add_node();
// 计算每个顶点与source和sink连接的权值，即第一个能量项
double sourceEdge, sinkEdge;
// 对于不确定的点，计算第一个能量项
if (mask.at<uchar>(i, j) == GC_PR_BGD || mask.at<uchar>(i, j) == GC_PR_FGD)
{
sourceEdge = -log(bgdGMM.GetProbability(img.at<Vec3b>(i, j)));
sinkEdge = -log(fgdGMM.GetProbability(img.at<Vec3b>(i, j)));
}
// 对于确定的背景点直接赋值
else if (mask.at<uchar>(i, j) == GC_BGD)
{
sourceEdge = 0;
sinkEdge = lambda;
}
else
{
sourceEdge = lambda;
sinkEdge = 0;
}
// 设置该node与source和sink点的权值
graph.add_tweights(nodeIndex, sourceEdge, sinkEdge);
// 左侧顶点的edge
if (j > 0)
{
double cap = leftWeight.at<double>(i, j);
graph.add_edge(nodeIndex, nodeIndex - 1, cap, cap);
}
// 上方顶点的edge
if (i > 0)
{
double cap = upWeight.at<double>(i, j);
graph.add_edge(nodeIndex, nodeIndex - img.cols, cap, cap);
}
// 左上方顶点的edge
if (j > 0 && i > 0)
{
double cap = upLeftWeight.at<double>(i, j);
graph.add_edge(nodeIndex, nodeIndex - img.cols - 1, cap, cap);
}
// 右上方顶点的edge
if (i > 0 && j < img.cols - 1)
{
double cap = upRightWeight.at<double>(i, j);
graph.add_edge(nodeIndex, nodeIndex - img.cols + 1, cap, cap);
}
}
}
}
// 使用最小割最大流算法进行分割估计
void GrabCut2D::EstimateSegmentation(Graph<double, double, double>& graph, Mat & mask)
{
graph.maxflow();
for (int i = 0; i < mask.rows; i++)
{
for (int j = 0; j < mask.cols; j++)
{
if (mask.at<uchar>(i, j) == GC_PR_BGD || mask.at<uchar>(i, j) == GC_PR_FGD)
{
if (graph.what_segment(i * mask.cols + j))
{
// 当返回1时为SINK
mask.at<uchar>(i, j) = GC_PR_BGD;
}
else
{
// 当返回0时为SOURCE
mask.at<uchar>(i, j) = GC_PR_FGD;
}
}
}
}
}
void GrabCut2D::GrabCut( cv::InputArray _img, cv::InputOutputArray _mask, cv::Rect rect, cv::InputOutputArray _bgdModel,cv::InputOutputArray _fgdModel, int iterCount, int mode )
{
std::cout<<"Execute GrabCut Function: Please finish the code here!"<<std::endl;
GMM bgdGmm(_bgdModel.getMatRef());
GMM fgdGmm(_fgdModel.getMatRef());
CheckRectAndMask(_mask.getMatRef(), rect, mode, _img.size());
InitGMMs(_img.getMat(), _mask.getMatRef(), bgdGmm, fgdGmm);
const double gamma = 50;
const double lambda = 9 * gamma;
const double beta = CalBeta(_img.getMat());
Mat leftWeight, upLeftWeight, upWeight, upRightWeight;
CalAllWeights(_img.getMat(), leftWeight, upLeftWeight, upWeight, upRightWeight, beta, gamma);
Mat compIdxs(_img.getMat().size(), CV_32SC1);
for (int i = 0; i < iterCount; i++)
{
int vertexCnt = _img.getMat().rows * _img.getMat().cols;
int edgeCnt = 2 * (4 * vertexCnt - 3 * (_img.getMat().cols + _img.getMat().rows) + 2);
Graph<double, double, double> graph(vertexCnt, edgeCnt);
AssignGMMsComponents(_img.getMat(), _mask.getMat(), bgdGmm, fgdGmm, compIdxs);
LearnGMMs(_img.getMat(), _mask.getMat(), compIdxs, bgdGmm, fgdGmm);
ConstructGraph(_img.getMat(), _mask.getMatRef(), bgdGmm, fgdGmm, lambda, leftWeight, upLeftWeight,
upWeight, upRightWeight, graph);
EstimateSegmentation(graph, _mask.getMatRef());
// cout << _mask.getMatRef() << endl;
}
//一.参数解释：
//输入：
//cv::InputArray _img, :输入的color图像(类型-cv:Mat)
//cv::Rect rect :在图像上画的矩形框（类型-cv:Rect)
//int iterCount : :每次分割的迭代次数（类型-int)
//中间变量
//cv::InputOutputArray _bgdModel ： 背景模型（推荐GMM)（类型-13*n（组件个数）个double类型的自定义数据结构，可以为cv:Mat，或者Vector/List/数组等）
//cv::InputOutputArray _fgdModel : 前景模型（推荐GMM) （类型-13*n（组件个数）个double类型的自定义数据结构，可以为cv:Mat，或者Vector/List/数组等）
//输出:
//cv::InputOutputArray _mask : 输出的分割结果 (类型： cv::Mat)
//二. 伪代码流程：
//1.Load Input Image: 加载输入颜色图像;
//2.Init Mask: 用矩形框初始化Mask的Label值（确定背景：0， 确定前景：1，可能背景：2，可能前景：3）,矩形框以外设置为确定背景，矩形框以内设置为可能前景;
//3.Init GMM: 定义并初始化GMM(其他模型完成分割也可得到基本分数，GMM完成会加分）
//4.Sample Points:前背景颜色采样并进行聚类（建议用kmeans，其他聚类方法也可)
//5.Learn GMM(根据聚类的样本更新每个GMM组件中的均值、协方差等参数）
//4.Construct Graph（计算t-weight(数据项）和n-weight（平滑项））
//7.Estimate Segmentation(调用maxFlow库进行分割)
//8.Save Result输入结果（将结果mask输出，将mask中前景区域对应的彩色图像保存和显示在交互界面中）
}
#pragma once
#include<opencv2\opencv.hpp>
#include <iostream>
#include "GMM.h"
#include "graph.h"
#include "block.h"
using namespace cv;
using namespace std;
enum
{
GC_WITH_RECT = 0,
GC_WITH_MASK = 1,
GC_CUT = 2
};
class GrabCut2D
{
public:
void GrabCut( cv::InputArray _img, cv::InputOutputArray _mask, cv::Rect rect,
cv::InputOutputArray _bgdModel,cv::InputOutputArray _fgdModel,
int iterCount, int mode );
void CheckRectAndMask(Mat &_mask, Rect rect, int mode, Size imgSize);
void InitMaskWithRect(Mat &_mask, Rect rect, Size imgSize);
void InitGMMs(Mat& img, Mat& mask, GMM& bgdGMM, GMM& fgdGMM);
double CalBeta(const Mat& img);
void CalAllWeights(const Mat& img, Mat& leftWeight, Mat& upLeftWeight, Mat& upWeight, Mat& upRightWeight,
double beta, double gamma);
void AssignGMMsComponents(const Mat& img, const Mat& mask, GMM& bgdGMM, GMM& fgdGMM,
Mat& compIdxs);
void LearnGMMs(const Mat& img, const Mat& mask, const Mat& compIdxs, GMM& bgdGMM, GMM& fgdGMM);
void ConstructGraph(const Mat& img, const Mat& mask, GMM& bgdGMM, GMM& fgdGMM,
double lambda, const Mat& leftWeight, const Mat& upLeftWeight, const Mat& upWeight,
const Mat upRightWeight, Graph<double, double, double>& graph);
void EstimateSegmentation(Graph<double, double, double>& graph, Mat& mask);
~GrabCut2D(void);
};
#pragma once
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "GrabCut.h"
#include <iostream>
#include "BorderMatting.h"
using namespace std;
using namespace cv;
const Scalar BLUE = Scalar(255,0,0); // Background
const Scalar GREEN = Scalar(0,255,0);//Foreground
const Scalar LIGHTBLUE = Scalar(255,255,160);//ProbBackground
const Scalar PINK = Scalar(230,130,255); //ProbBackground
const Scalar RED = Scalar(0,0,255);//color of Rectangle
const int BGD_KEY = CV_EVENT_FLAG_CTRLKEY;// When press "CTRL" key,the value of flags return.
const int FGD_KEY = CV_EVENT_FLAG_SHIFTKEY;// When press "SHIFT" key, the value of flags return.
//Copy the value of comMask to binMask
static void getBinMask( const Mat& comMask, Mat& binMask )
{
if( comMask.empty() || comMask.type()!=CV_8UC1 )
CV_Error( CV_StsBadArg, "comMask is empty or has incorrect type (not CV_8UC1)" );
if( binMask.empty() || binMask.rows!=comMask.rows || binMask.cols!=comMask.cols )
binMask.create( comMask.size(), CV_8UC1 );
binMask = comMask & 1;
}
class GCApplication
{
public:
enum{ NOT_SET = 0, IN_PROCESS = 1, SET = 2 };
static const int radius = 2;
static const int thickness = -1;
void reset();
void setImageAndWinName( const Mat& _image, const string& _winName );
void showImage() const;
void mouseClick( int event, int x, int y, int flags, void* param );
int nextIter();
int getIterCount() const { return iterCount; }
void borderMatting();
private:
void setRectInMask();
void setLblsInMask( int flags, Point p, bool isPr );
const string* winName;
const Mat* image;
Mat mask, alphaMask;
Mat bgdModel, fgdModel;
uchar rectState, lblsState, prLblsState;
bool isInitialized;
Rect rect;
vector<Point> fgdPxls, bgdPxls, prFgdPxls, prBgdPxls;
int iterCount;
GrabCut2D gc;
BorderMatting bm;
};
#include "GCApplication.h"
//Set value for the class
void GCApplication::reset()
{
if( !mask.empty() )
mask.setTo(Scalar::all(GC_BGD));
bgdPxls.clear(); fgdPxls.clear();
prBgdPxls.clear(); prFgdPxls.clear();
isInitialized = false;
rectState = NOT_SET;
lblsState = NOT_SET;
prLblsState = NOT_SET;
iterCount = 0;
}
//Set image and window name
void GCApplication::setImageAndWinName( const Mat& _image, const string& _winName )
{
if( _image.empty() || _winName.empty() )
return;
image = &_image;
winName = &_winName;
mask.create( image->size(), CV_8UC1);
reset();
}
//Show the result image
void GCApplication::showImage() const
{
if( image->empty() || winName->empty() )
return;
Mat res;
Mat binMask;
if( !isInitialized )
image->copyTo( res );
else
{
getBinMask( mask, binMask );
image->copyTo( res, binMask ); //show the GrabCuted image
}
vector<Point>::const_iterator it;
//Using four different colors show the point which have been selected
for( it = bgdPxls.begin(); it != bgdPxls.end(); ++it )
circle( res, *it, radius, BLUE, thickness );
for( it = fgdPxls.begin(); it != fgdPxls.end(); ++it )
circle( res, *it, radius, GREEN, thickness );
for( it = prBgdPxls.begin(); it != prBgdPxls.end(); ++it )
circle( res, *it, radius, LIGHTBLUE, thickness );
for( it = prFgdPxls.begin(); it != prFgdPxls.end(); ++it )
circle( res, *it, radius, PINK, thickness );
imwrite("result.jpg", res);
//Draw the rectangle
if( rectState == IN_PROCESS || rectState == SET )
rectangle( res, Point( rect.x, rect.y ), Point(rect.x + rect.width, rect.y + rect.height ), RED, 2);
imshow( *winName, res );
}
//Using rect initialize the pixel
void GCApplication::setRectInMask()
{
assert( !mask.empty() );
mask.setTo( GC_BGD ); //GC_BGD == 0
rect.x = max(0, rect.x);
rect.y = max(0, rect.y);
rect.width = min(rect.width, image->cols-rect.x);
rect.height = min(rect.height, image->rows-rect.y);
(mask(rect)).setTo( Scalar(GC_PR_FGD) ); //GC_PR_FGD == 3
}
void GCApplication::setLblsInMask( int flags, Point p, bool isPr )
{
vector<Point> *bpxls, *fpxls;
uchar bvalue, fvalue;
if( !isPr ) //Points which are sure being FGD or BGD
{
bpxls = &bgdPxls;
fpxls = &fgdPxls;
bvalue = GC_BGD; //0
fvalue = GC_FGD; //1
}
else //Probably FGD or Probably BGD
{
bpxls = &prBgdPxls;
fpxls = &prFgdPxls;
bvalue = GC_PR_BGD; //2
fvalue = GC_PR_FGD; //3
}
if( flags & BGD_KEY )
{
bpxls->push_back(p);
circle( mask, p, radius, bvalue, thickness ); //Set point value = 2
}
if( flags & FGD_KEY )
{
fpxls->push_back(p);
circle( mask, p, radius, fvalue, thickness ); //Set point value = 3
}
}
//Mouse Click Function: flags work with CV_EVENT_FLAG
void GCApplication::mouseClick( int event, int x, int y, int flags, void* )
{
switch( event )
{
case CV_EVENT_LBUTTONDOWN: // Set rect or GC_BGD(GC_FGD) labels
{
bool isb = (flags & BGD_KEY) != 0,
isf = (flags & FGD_KEY) != 0;
if( rectState == NOT_SET && !isb && !isf )//Only LEFT_KEY pressed
{
rectState = IN_PROCESS; //Be drawing the rectangle
rect = Rect( x, y, 1, 1 );
}
if ( (isb || isf) && rectState == SET ) //Set the BGD/FGD(labels).after press the "ALT" key or "SHIFT" key,and have finish drawing the rectangle
lblsState = IN_PROCESS;
}
break;
case CV_EVENT_RBUTTONDOWN: // Set GC_PR_BGD(GC_PR_FGD) labels
{
bool isb = (flags & BGD_KEY) != 0,
isf = (flags & FGD_KEY) != 0;
if ( (isb || isf) && rectState == SET ) //Set the probably FGD/BGD labels
prLblsState = IN_PROCESS;
}
break;
case CV_EVENT_LBUTTONUP:
if( rectState == IN_PROCESS )
{
rect = Rect( Point(rect.x, rect.y), Point(x,y) ); //After draw the rectangle
rectState = SET;
setRectInMask();
assert( bgdPxls.empty() && fgdPxls.empty() && prBgdPxls.empty() && prFgdPxls.empty() );
showImage();
}
if( lblsState == IN_PROCESS )
{
setLblsInMask(flags, Point(x,y), false); // Draw the FGD points
lblsState = SET;
showImage();
}
break;
case CV_EVENT_RBUTTONUP:
if( prLblsState == IN_PROCESS )
{
setLblsInMask(flags, Point(x,y), true); //Draw the BGD points
prLblsState = SET;
showImage();
}
break;
case CV_EVENT_MOUSEMOVE:
if( rectState == IN_PROCESS )
{
rect = Rect( Point(rect.x, rect.y), Point(x,y) );
assert( bgdPxls.empty() && fgdPxls.empty() && prBgdPxls.empty() && prFgdPxls.empty() );
showImage(); //Continue showing image
}
else if( lblsState == IN_PROCESS )
{
setLblsInMask(flags, Point(x,y), false);
showImage();
}
else if( prLblsState == IN_PROCESS )
{
setLblsInMask(flags, Point(x,y), true);
showImage();
}
break;
}
}
void GCApplication::borderMatting() {
bm.BorderMattingInterface(*image, mask, alphaMask);
}
//Execute GrabCut algorithm，and return the iter count.
int GCApplication::nextIter()
{
if( isInitialized )
gc.GrabCut(*image, mask, rect, bgdModel, fgdModel,1,GC_CUT);
else
{
if( rectState != SET )
return iterCount;
if( lblsState == SET || prLblsState == SET )
gc.GrabCut( *image, mask, rect, bgdModel, fgdModel, 1, GC_WITH_MASK );
else
gc.GrabCut(*image, mask, rect, bgdModel, fgdModel,1,GC_WITH_RECT);
isInitialized = true;
}
iterCount++;
bgdPxls.clear(); fgdPxls.clear();
prBgdPxls.clear(); prFgdPxls.clear();
return iterCount;
}
#include <iostream>
#include "GCApplication.h"
static void help()
{
std::cout << "\nThis program demonstrates GrabCut segmentation -- select an object in a region\n"
"and then grabcut will attempt to segment it out.\n"
"Call:\n"
"./grabcut <image_name>\n"
"\nSelect a rectangular area around the object you want to segment\n" <<
"\nHot keys: \n"
"\tESC - quit the program\n"
"\tr - restore the original image\n"
"\tn - next iteration\n"
"\n"
"\tleft mouse button - set rectangle\n"
"\n"
"\tCTRL+left mouse button - set GC_BGD pixels\n"
"\tSHIFT+left mouse button - set CG_FGD pixels\n"
"\n"
"\tCTRL+right mouse button - set GC_PR_BGD pixels\n"
"\tSHIFT+right mouse button - set CG_PR_FGD pixels\n" << endl;
}
GCApplication gcapp;
static void on_mouse( int event, int x, int y, int flags, void* param )
{
gcapp.mouseClick( event, x, y, flags, param );
}
int main()
{
string filename = "img.jpg";
Mat image = imread( filename, 1 );
//resize(image, image, Size(image.cols / 4, image.rows / 4));
if( image.empty() )
{
cout << "\n , couldn't read image filename " << filename << endl;
return 1;
}
help();
const string winName = "image";
cvNamedWindow( winName.c_str(), CV_WINDOW_AUTOSIZE );
cvSetMouseCallback( winName.c_str(), on_mouse, 0 );
gcapp.setImageAndWinName( image, winName );
gcapp.showImage();
for(;;)
{
int c = cvWaitKey(0);
switch( (char) c )
{
case '\x1b':
cout << "Exiting ..." << endl;
goto exit_main;
case 'r':
cout << endl;
gcapp.reset();
gcapp.showImage();
break;
case 'b':
gcapp.borderMatting();
cout << "done" << endl;
break;
case 'n':
int iterCount = gcapp.getIterCount();
cout << "<" << iterCount << "... ";
int newIterCount = gcapp.nextIter();
if( newIterCount > iterCount )
{
gcapp.showImage();
cout << iterCount << ">" << endl;
}
else
cout << "rect must be determined>" << endl;
break;
}
}
exit_main:
cvDestroyWindow( winName.c_str() );
return 0;
return 0;
}
/* graph.h */
/*
This software library implements the maxflow algorithm
described in
"An Experimental Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision."
Yuri Boykov and Vladimir Kolmogorov.
In IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI),
September 2004
This algorithm was developed by Yuri Boykov and Vladimir Kolmogorov
at Siemens Corporate Research. To make it available for public use,
it was later reimplemented by Vladimir Kolmogorov based on open publications.
If you use this software for research purposes, you should cite
the aforementioned paper in any resulting publication.
----------------------------------------------------------------------
REUSING TREES:
Starting with version 3.0, there is a also an option of reusing search
trees from one maxflow computation to the next, as described in
"Efficiently Solving Dynamic Markov Random Fields Using Graph Cuts."
Pushmeet Kohli and Philip H.S. Torr
International Conference on Computer Vision (ICCV), 2005
If you use this option, you should cite
the aforementioned paper in any resulting publication.
*/
/*
For description, license, example usage see README.TXT.
*/
#ifndef __GRAPH_H__
#define __GRAPH_H__
#include <string.h>
#include "block.h"
#include <assert.h>
// NOTE: in UNIX you need to use -DNDEBUG preprocessor option to supress assert's!!!
// captype: type of edge capacities (excluding t-links)
// tcaptype: type of t-links (edges between nodes and terminals)
// flowtype: type of total flow
//
// Current instantiations are in instances.inc
template <typename captype, typename tcaptype, typename flowtype> class Graph
{
public:
typedef enum
{
SOURCE = 0,
SINK = 1
} termtype; // terminals
typedef int node_id;
/
// BASIC INTERFACE FUNCTIONS //
// (should be enough for most applications) //
/
// Constructor.
// The first argument gives an estimate of the maximum number of nodes that can be added
// to the graph, and the second argument is an estimate of the maximum number of edges.
// The last (optional) argument is the pointer to the function which will be called
// if an error occurs; an error message is passed to this function.
// If this argument is omitted, exit(1) will be called.
//
// IMPORTANT: It is possible to add more nodes to the graph than node_num_max
// (and node_num_max can be zero). However, if the count is exceeded, then
// the internal memory is reallocated (increased by 50%) which is expensive.
// Also, temporarily the amount of allocated memory would be more than twice than needed.
// Similarly for edges.
// If you wish to avoid this overhead, you can download version 2.2, where nodes and edges are stored in blocks.
Graph(int node_num_max, int edge_num_max, void (*err_function)(char *) = NULL);
// Destructor
~Graph();
// Adds node(s) to the graph. By default, one node is added (num=1); then first call returns 0, second call returns 1, and so on.
// If num>1, then several nodes are added, and node_id of the first one is returned.
// IMPORTANT: see note about the constructor
node_id add_node(int num = 1);
// Adds a bidirectional edge between 'i' and 'j' with the weights 'cap' and 'rev_cap'.
// IMPORTANT: see note about the constructor
void add_edge(node_id i, node_id j, captype cap, captype rev_cap);
// Adds new edges 'SOURCE->i' and 'i->SINK' with corresponding weights.
// Can be called multiple times for each node.
// Weights can be negative.
// NOTE: the number of such edges is not counted in edge_num_max.
// No internal memory is allocated by this call.
void add_tweights(node_id i, tcaptype cap_source, tcaptype cap_sink);
// Computes the maxflow. Can be called several times.
// FOR DESCRIPTION OF reuse_trees, SEE mark_node().
// FOR DESCRIPTION OF changed_list, SEE remove_from_changed_list().
flowtype maxflow(bool reuse_trees = false, Block<node_id>* changed_list = NULL);
// After the maxflow is computed, this function returns to which
// segment the node 'i' belongs (Graph<captype,tcaptype,flowtype>::SOURCE or Graph<captype,tcaptype,flowtype>::SINK).
//
// Occasionally there may be several minimum cuts. If a node can be assigned
// to both the source and the sink, then default_segm is returned.
termtype what_segment(node_id i, termtype default_segm = SOURCE);
//
// ADVANCED INTERFACE FUNCTIONS //
// (provide access to the graph) //
//
private:
struct node;
struct arc;
public:
// 1. Reallocating graph. //
// Removes all nodes and edges.
// After that functions add_node() and add_edge() must be called again.
//
// Advantage compared to deleting Graph and allocating it again:
// no calls to delete/new (which could be quite slow).
//
// If the graph structure stays the same, then an alternative
// is to go through all nodes/edges and set new residual capacities
// (see functions below).
void reset();
// 2. Functions for getting pointers to arcs and for reading graph structure. //
// NOTE: adding new arcs may invalidate these pointers (if reallocation //
// happens). So it's best not to add arcs while reading graph structure. //
// The following two functions return arcs in the same order that they
// were added to the graph. NOTE: for each call add_edge(i,j,cap,cap_rev)
// the first arc returned will be i->j, and the second j->i.
// If there are no more arcs, then the function can still be called, but
// the returned arc_id is undetermined.
typedef arc* arc_id;
arc_id get_first_arc();
arc_id get_next_arc(arc_id a);
// other functions for reading graph structure
int get_node_num() { return node_num; }
int get_arc_num() { return (int)(arc_last - arcs); }
void get_arc_ends(arc_id a, node_id& i, node_id& j); // returns i,j to that a = i->j
///
// 3. Functions for reading residual capacities. //
///
// returns residual capacity of SOURCE->i minus residual capacity of i->SINK
tcaptype get_trcap(node_id i);
// returns residual capacity of arc a
captype get_rcap(arc* a);
/
// 4. Functions for setting residual capacities. //
// NOTE: If these functions are used, the value of the flow //
// returned by maxflow() will not be valid! //
/
void set_trcap(node_id i, tcaptype trcap);
void set_rcap(arc* a, captype rcap);
// 5. Functions related to reusing trees & list of changed nodes. //
// If flag reuse_trees is true while calling maxflow(), then search trees
// are reused from previous maxflow computation.
// In this case before calling maxflow() the user must
// specify which parts of the graph have changed by calling mark_node():
// add_tweights(i),set_trcap(i) => call mark_node(i)
// add_edge(i,j),set_rcap(a) => call mark_node(i); mark_node(j)
//
// This option makes sense only if a small part of the graph is changed.
// The initialization procedure goes only through marked nodes then.
//
// mark_node(i) can either be called before or after graph modification.
// Can be called more than once per node, but calls after the first one
// do not have any effect.
//
// NOTE:
// - This option cannot be used in the first call to maxflow().
// - It is not necessary to call mark_node() if the change is ``not essential'',
// i.e. sign(trcap) is preserved for a node and zero/nonzero status is preserved for an arc.
// - To check that you marked all necessary nodes, you can call maxflow(false) after calling maxflow(true).
// If everything is correct, the two calls must return the same value of flow. (Useful for debugging).
void mark_node(node_id i);
// If changed_list is not NULL while calling maxflow(), then the algorithm
// keeps a list of nodes which could potentially have changed their segmentation label.
// Nodes which are not in the list are guaranteed to keep their old segmentation label (SOURCE or SINK).
// Example usage:
//
// typedef Graph<int,int,int> G;
// G* g = new Graph(nodeNum, edgeNum);
// Block<G::node_id>* changed_list = new Block<G::node_id>(128);
//
// ... // add nodes and edges
//
// g->maxflow(); // first call should be without arguments
// for (int iter=0; iter<10; iter++)
// {
// ... // change graph, call mark_node() accordingly
//
// g->maxflow(true, changed_list);
// G::node_id* ptr;
// for (ptr=changed_list->ScanFirst(); ptr; ptr=changed_list->ScanNext())
// {
// G::node_id i = *ptr; assert(i>=0 && i<nodeNum);
// g->remove_from_changed_list(i);
// // do something with node i...
// if (g->what_segment(i) == G::SOURCE) { ... }
// }
// changed_list->Reset();
// }
// delete changed_list;
//
// NOTE:
// - If changed_list option is used, then reuse_trees must be used as well.
// - In the example above, the user may omit calls g->remove_from_changed_list(i) and changed_list->Reset() in a given iteration.
// Then during the next call to maxflow(true, &changed_list) new nodes will be added to changed_list.
// - If the next call to maxflow() does not use option reuse_trees, then calling remove_from_changed_list()
// is not necessary. ("changed_list->Reset()" or "delete changed_list" should still be called, though).
void remove_from_changed_list(node_id i)
{
assert(i>=0 && i<node_num && nodes[i].is_in_changed_list);
nodes[i].is_in_changed_list = 0;
}
/
/
/
private:
// internal variables and functions
struct node
{
arc *first; // first outcoming arc
arc *parent; // node's parent
node *next; // pointer to the next active node
// (or to itself if it is the last node in the list)
int TS; // timestamp showing when DIST was computed
int DIST; // distance to the terminal
int is_sink : 1; // flag showing whether the node is in the source or in the sink tree (if parent!=NULL)
int is_marked : 1; // set by mark_node()
int is_in_changed_list : 1; // set by maxflow if
tcaptype tr_cap; // if tr_cap > 0 then tr_cap is residual capacity of the arc SOURCE->node
// otherwise -tr_cap is residual capacity of the arc node->SINK
};
struct arc
{
node *head; // node the arc points to
arc *next; // next arc with the same originating node
arc *sister; // reverse arc
captype r_cap; // residual capacity
};
struct nodeptr
{
node *ptr;
nodeptr *next;
};
static const int NODEPTR_BLOCK_SIZE = 128;
node *nodes, *node_last, *node_max; // node_last = nodes+node_num, node_max = nodes+node_num_max;
arc *arcs, *arc_last, *arc_max; // arc_last = arcs+2*edge_num, arc_max = arcs+2*edge_num_max;
int node_num;
DBlock<nodeptr> *nodeptr_block;
void (*error_function)(char *); // this function is called if a error occurs,
// with a corresponding error message
// (or exit(1) is called if it's NULL)
flowtype flow; // total flow
// reusing trees & list of changed pixels
int maxflow_iteration; // counter
Block<node_id> *changed_list;
/
node *queue_first[2], *queue_last[2]; // list of active nodes
nodeptr *orphan_first, *orphan_last; // list of pointers to orphans
int TIME; // monotonically increasing global counter
/
void reallocate_nodes(int num); // num is the number of new nodes
void reallocate_arcs();
// functions for processing active list
void set_active(node *i);
node *next_active();
// functions for processing orphans list
void set_orphan_front(node* i); // add to the beginning of the list
void set_orphan_rear(node* i); // add to the end of the list
void add_to_changed_list(node* i);
void maxflow_init(); // called if reuse_trees == false
void maxflow_reuse_trees_init(); // called if reuse_trees == true
void augment(arc *middle_arc);
void process_source_orphan(node *i);
void process_sink_orphan(node *i);
void test_consistency(node* current_node=NULL); // debug function
};
///
// Implementation - inline functions //
///
template <typename captype, typename tcaptype, typename flowtype>
inline typename Graph<captype,tcaptype,flowtype>::node_id Graph<captype,tcaptype,flowtype>::add_node(int num)
{
assert(num > 0);
if (node_last + num > node_max) reallocate_nodes(num);
if (num == 1)
{
node_last -> first = NULL;
node_last -> tr_cap = 0;
node_last -> is_marked = 0;
node_last -> is_in_changed_list = 0;
node_last ++;
return node_num ++;
}
else
{
memset(node_last, 0, num*sizeof(node));
node_id i = node_num;
node_num += num;
node_last += num;
return i;
}
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::add_tweights(node_id i, tcaptype cap_source, tcaptype cap_sink)
{
assert(i >= 0 && i < node_num);
tcaptype delta = nodes[i].tr_cap;
if (delta > 0) cap_source += delta;
else cap_sink -= delta;
flow += (cap_source < cap_sink) ? cap_source : cap_sink;
nodes[i].tr_cap = cap_source - cap_sink;
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::add_edge(node_id _i, node_id _j, captype cap, captype rev_cap)
{
assert(_i >= 0 && _i < node_num);
assert(_j >= 0 && _j < node_num);
assert(_i != _j);
assert(cap >= 0);
assert(rev_cap >= 0);
if (arc_last == arc_max) reallocate_arcs();
arc *a = arc_last ++;
arc *a_rev = arc_last ++;
node* i = nodes + _i;
node* j = nodes + _j;
a -> sister = a_rev;
a_rev -> sister = a;
a -> next = i -> first;
i -> first = a;
a_rev -> next = j -> first;
j -> first = a_rev;
a -> head = j;
a_rev -> head = i;
a -> r_cap = cap;
a_rev -> r_cap = rev_cap;
}
template <typename captype, typename tcaptype, typename flowtype>
inline typename Graph<captype,tcaptype,flowtype>::arc* Graph<captype,tcaptype,flowtype>::get_first_arc()
{
return arcs;
}
template <typename captype, typename tcaptype, typename flowtype>
inline typename Graph<captype,tcaptype,flowtype>::arc* Graph<captype,tcaptype,flowtype>::get_next_arc(arc* a)
{
return a + 1;
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::get_arc_ends(arc* a, node_id& i, node_id& j)
{
assert(a >= arcs && a < arc_last);
i = (node_id) (a->sister->head - nodes);
j = (node_id) (a->head - nodes);
}
template <typename captype, typename tcaptype, typename flowtype>
inline tcaptype Graph<captype,tcaptype,flowtype>::get_trcap(node_id i)
{
assert(i>=0 && i<node_num);
return nodes[i].tr_cap;
}
template <typename captype, typename tcaptype, typename flowtype>
inline captype Graph<captype,tcaptype,flowtype>::get_rcap(arc* a)
{
assert(a >= arcs && a < arc_last);
return a->r_cap;
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::set_trcap(node_id i, tcaptype trcap)
{
assert(i>=0 && i<node_num);
nodes[i].tr_cap = trcap;
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::set_rcap(arc* a, captype rcap)
{
assert(a >= arcs && a < arc_last);
a->r_cap = rcap;
}
template <typename captype, typename tcaptype, typename flowtype>
inline typename Graph<captype,tcaptype,flowtype>::termtype Graph<captype,tcaptype,flowtype>::what_segment(node_id i, termtype default_segm)
{
if (nodes[i].parent)
{
return (nodes[i].is_sink) ? SINK : SOURCE;
}
else
{
return default_segm;
}
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::mark_node(node_id _i)
{
node* i = nodes + _i;
if (!i->next)
{
/* it's not in the list yet */
if (queue_last[1]) queue_last[1] -> next = i;
else queue_first[1] = i;
queue_last[1] = i;
i -> next = i;
}
i->is_marked = 1;
}
#endif
/* graph.cpp */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "graph.h"
template <typename captype, typename tcaptype, typename flowtype>
Graph<captype, tcaptype, flowtype>::Graph(int node_num_max, int edge_num_max, void (*err_function)(char *))
: node_num(0),
nodeptr_block(NULL),
error_function(err_function)
{
if (node_num_max < 16) node_num_max = 16;
if (edge_num_max < 16) edge_num_max = 16;
nodes = (node*) malloc(node_num_max*sizeof(node));
arcs = (arc*) malloc(2*edge_num_max*sizeof(arc));
if (!nodes || !arcs) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
node_last = nodes;
node_max = nodes + node_num_max;
arc_last = arcs;
arc_max = arcs + 2*edge_num_max;
maxflow_iteration = 0;
flow = 0;
}
template <typename captype, typename tcaptype, typename flowtype>
Graph<captype,tcaptype,flowtype>::~Graph()
{
if (nodeptr_block)
{
delete nodeptr_block;
nodeptr_block = NULL;
}
free(nodes);
free(arcs);
}
template <typename captype, typename tcaptype, typename flowtype>
void Graph<captype,tcaptype,flowtype>::reset()
{
node_last = nodes;
arc_last = arcs;
node_num = 0;
if (nodeptr_block)
{
delete nodeptr_block;
nodeptr_block = NULL;
}
maxflow_iteration = 0;
flow = 0;
}
template <typename captype, typename tcaptype, typename flowtype>
void Graph<captype,tcaptype,flowtype>::reallocate_nodes(int num)
{
int node_num_max = (int)(node_max - nodes);
node* nodes_old = nodes;
node_num_max += node_num_max / 2;
if (node_num_max < node_num + num) node_num_max = node_num + num;
nodes = (node*) realloc(nodes_old, node_num_max*sizeof(node));
if (!nodes) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
node_last = nodes + node_num;
node_max = nodes + node_num_max;
if (nodes != nodes_old)
{
arc* a;
for (a=arcs; a<arc_last; a++)
{
a->head = (node*) ((char*)a->head + (((char*) nodes) - ((char*) nodes_old)));
}
}
}
template <typename captype, typename tcaptype, typename flowtype>
void Graph<captype,tcaptype,flowtype>::reallocate_arcs()
{
int arc_num_max = (int)(arc_max - arcs);
int arc_num = (int)(arc_last - arcs);
arc* arcs_old = arcs;
arc_num_max += arc_num_max / 2; if (arc_num_max & 1) arc_num_max ++;
arcs = (arc*) realloc(arcs_old, arc_num_max*sizeof(arc));
if (!arcs) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
arc_last = arcs + arc_num;
arc_max = arcs + arc_num_max;
if (arcs != arcs_old)
{
node* i;
arc* a;
for (i=nodes; i<node_last; i++)
{
if (i->first) i->first = (arc*) ((char*)i->first + (((char*) arcs) - ((char*) arcs_old)));
}
for (a=arcs; a<arc_last; a++)
{
if (a->next) a->next = (arc*) ((char*)a->next + (((char*) arcs) - ((char*) arcs_old)));
a->sister = (arc*) ((char*)a->sister + (((char*) arcs) - ((char*) arcs_old)));
}
}
}
#include "instances.inc"
/* maxflow.cpp */
#include <stdio.h>
#include "graph.h"
/*
special constants for node->parent
*/
#define TERMINAL ( (arc *) 1 ) /* to terminal */
#define ORPHAN ( (arc *) 2 ) /* orphan */
#define INFINITE_D ((int)(((unsigned)-1)/2)) /* infinite distance to the terminal */
/***********************************************************************/
/*
Functions for processing active list.
i->next points to the next node in the list
(or to i, if i is the last node in the list).
If i->next is NULL iff i is not in the list.
There are two queues. Active nodes are added
to the end of the second queue and read from
the front of the first queue. If the first queue
is empty, it is replaced by the second queue
(and the second queue becomes empty).
*/
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::set_active(node *i)
{
if (!i->next)
{
/* it's not in the list yet */
if (queue_last[1]) queue_last[1] -> next = i;
else queue_first[1] = i;
queue_last[1] = i;
i -> next = i;
}
}
/*
Returns the next active node.
If it is connected to the sink, it stays in the list,
otherwise it is removed from the list
*/
template <typename captype, typename tcaptype, typename flowtype>
inline typename Graph<captype,tcaptype,flowtype>::node* Graph<captype,tcaptype,flowtype>::next_active()
{
node *i;
while ( 1 )
{
if (!(i=queue_first[0]))
{
queue_first[0] = i = queue_first[1];
queue_last[0] = queue_last[1];
queue_first[1] = NULL;
queue_last[1] = NULL;
if (!i) return NULL;
}
/* remove it from the active list */
if (i->next == i) queue_first[0] = queue_last[0] = NULL;
else queue_first[0] = i -> next;
i -> next = NULL;
/* a node in the list is active iff it has a parent */
if (i->parent) return i;
}
}
/***********************************************************************/
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::set_orphan_front(node *i)
{
nodeptr *np;
i -> parent = ORPHAN;
np = nodeptr_block -> New();
np -> ptr = i;
np -> next = orphan_first;
orphan_first = np;
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::set_orphan_rear(node *i)
{
nodeptr *np;
i -> parent = ORPHAN;
np = nodeptr_block -> New();
np -> ptr = i;
if (orphan_last) orphan_last -> next = np;
else orphan_first = np;
orphan_last = np;
np -> next = NULL;
}
/***********************************************************************/
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::add_to_changed_list(node *i)
{
if (changed_list && !i->is_in_changed_list)
{
node_id* ptr = changed_list->New();
*ptr = (node_id)(i - nodes);
i->is_in_changed_list = true;
}
}
/***********************************************************************/
template <typename captype, typename tcaptype, typename flowtype>
void Graph<captype,tcaptype,flowtype>::maxflow_init()
{
node *i;
queue_first[0] = queue_last[0] = NULL;
queue_first[1] = queue_last[1] = NULL;
orphan_first = NULL;
TIME = 0;
for (i=nodes; i<node_last; i++)
{
i -> next = NULL;
i -> is_marked = 0;
i -> is_in_changed_list = 0;
i -> TS = TIME;
if (i->tr_cap > 0)
{
/* i is connected to the source */
i -> is_sink = 0;
i -> parent = TERMINAL;
set_active(i);
i -> DIST = 1;
}
else if (i->tr_cap < 0)
{
/* i is connected to the sink */
i -> is_sink = 1;
i -> parent = TERMINAL;
set_active(i);
i -> DIST = 1;
}
else
{
i -> parent = NULL;
}
}
}
template <typename captype, typename tcaptype, typename flowtype>
void Graph<captype,tcaptype,flowtype>::maxflow_reuse_trees_init()
{
node* i;
node* j;
node* queue = queue_first[1];
arc* a;
nodeptr* np;
queue_first[0] = queue_last[0] = NULL;
queue_first[1] = queue_last[1] = NULL;
orphan_first = orphan_last = NULL;
TIME ++;
while ((i=queue))
{
queue = i->next;
if (queue == i) queue = NULL;
i->next = NULL;
i->is_marked = 0;
set_active(i);
if (i->tr_cap == 0)
{
if (i->parent) set_orphan_rear(i);
continue;
}
if (i->tr_cap > 0)
{
if (!i->parent || i->is_sink)
{
i->is_sink = 0;
for (a=i->first; a; a=a->next)
{
j = a->head;
if (!j->is_marked)
{
if (j->parent == a->sister) set_orphan_rear(j);
if (j->parent && j->is_sink && a->r_cap > 0) set_active(j);
}
}
add_to_changed_list(i);
}
}
else
{
if (!i->parent || !i->is_sink)
{
i->is_sink = 1;
for (a=i->first; a; a=a->next)
{
j = a->head;
if (!j->is_marked)
{
if (j->parent == a->sister) set_orphan_rear(j);
if (j->parent && !j->is_sink && a->sister->r_cap > 0) set_active(j);
}
}
add_to_changed_list(i);
}
}
i->parent = TERMINAL;
i -> TS = TIME;
i -> DIST = 1;
}
//test_consistency();
/* adoption */
while ((np=orphan_first))
{
orphan_first = np -> next;
i = np -> ptr;
nodeptr_block -> Delete(np);
if (!orphan_first) orphan_last = NULL;
if (i->is_sink) process_sink_orphan(i);
else process_source_orphan(i);
}
/* adoption end */
//test_consistency();
}
template <typename captype, typename tcaptype, typename flowtype>
void Graph<captype,tcaptype,flowtype>::augment(arc *middle_arc)
{
node *i;
arc *a;
tcaptype bottleneck;
/* 1. Finding bottleneck capacity */
/* 1a - the source tree */
bottleneck = middle_arc -> r_cap;
for (i=middle_arc->sister->head; ; i=a->head)
{
a = i -> parent;
if (a == TERMINAL) break;
if (bottleneck > a->sister->r_cap) bottleneck = a -> sister -> r_cap;
}
if (bottleneck > i->tr_cap) bottleneck = i -> tr_cap;
/* 1b - the sink tree */
for (i=middle_arc->head; ; i=a->head)
{
a = i -> parent;
if (a == TERMINAL) break;
if (bottleneck > a->r_cap) bottleneck = a -> r_cap;
}
if (bottleneck > - i->tr_cap) bottleneck = - i -> tr_cap;
/* 2. Augmenting */
/* 2a - the source tree */
middle_arc -> sister -> r_cap += bottleneck;
middle_arc -> r_cap -= bottleneck;
for (i=middle_arc->sister->head; ; i=a->head)
{
a = i -> parent;
if (a == TERMINAL) break;
a -> r_cap += bottleneck;
a -> sister -> r_cap -= bottleneck;
if (!a->sister->r_cap)
{
set_orphan_front(i); // add i to the beginning of the adoption list
}
}
i -> tr_cap -= bottleneck;
if (!i->tr_cap)
{
set_orphan_front(i); // add i to the beginning of the adoption list
}
/* 2b - the sink tree */
for (i=middle_arc->head; ; i=a->head)
{
a = i -> parent;
if (a == TERMINAL) break;
a -> sister -> r_cap += bottleneck;
a -> r_cap -= bottleneck;
if (!a->r_cap)
{
set_orphan_front(i); // add i to the beginning of the adoption list
}
}
i -> tr_cap += bottleneck;
if (!i->tr_cap)
{
set_orphan_front(i); // add i to the beginning of the adoption list
}
flow += bottleneck;
}
/***********************************************************************/
template <typename captype, typename tcaptype, typename flowtype>
void Graph<captype,tcaptype,flowtype>::process_source_orphan(node *i)
{
node *j;
arc *a0, *a0_min = NULL, *a;
int d, d_min = INFINITE_D;
/* trying to find a new parent */
for (a0=i->first; a0; a0=a0->next)
if (a0->sister->r_cap)
{
j = a0 -> head;
if (!j->is_sink && (a=j->parent))
{
/* checking the origin of j */
d = 0;
while ( 1 )
{
if (j->TS == TIME)
{
d += j -> DIST;
break;
}
a = j -> parent;
d ++;
if (a==TERMINAL)
{
j -> TS = TIME;
j -> DIST = 1;
break;
}
if (a==ORPHAN) { d = INFINITE_D; break; }
j = a -> head;
}
if (d<INFINITE_D) /* j originates from the source - done */
{
if (d<d_min)
{
a0_min = a0;
d_min = d;
}
/* set marks along the path */
for (j=a0->head; j->TS!=TIME; j=j->parent->head)
{
j -> TS = TIME;
j -> DIST = d --;
}
}
}
}
if (i->parent = a0_min)
{
i -> TS = TIME;
i -> DIST = d_min + 1;
}
else
{
/* no parent is found */
add_to_changed_list(i);
/* process neighbors */
for (a0=i->first; a0; a0=a0->next)
{
j = a0 -> head;
if (!j->is_sink && (a=j->parent))
{
if (a0->sister->r_cap) set_active(j);
if (a!=TERMINAL && a!=ORPHAN && a->head==i)
{
set_orphan_rear(j); // add j to the end of the adoption list
}
}
}
}
}
template <typename captype, typename tcaptype, typename flowtype>
void Graph<captype,tcaptype,flowtype>::process_sink_orphan(node *i)
{
node *j;
arc *a0, *a0_min = NULL, *a;
int d, d_min = INFINITE_D;
/* trying to find a new parent */
for (a0=i->first; a0; a0=a0->next)
if (a0->r_cap)
{
j = a0 -> head;
if (j->is_sink && (a=j->parent))
{
/* checking the origin of j */
d = 0;
while ( 1 )
{
if (j->TS == TIME)
{
d += j -> DIST;
break;
}
a = j -> parent;
d ++;
if (a==TERMINAL)
{
j -> TS = TIME;
j -> DIST = 1;
break;
}
if (a==ORPHAN) { d = INFINITE_D; break; }
j = a -> head;
}
if (d<INFINITE_D) /* j originates from the sink - done */
{
if (d<d_min)
{
a0_min = a0;
d_min = d;
}
/* set marks along the path */
for (j=a0->head; j->TS!=TIME; j=j->parent->head)
{
j -> TS = TIME;
j -> DIST = d --;
}
}
}
}
if (i->parent = a0_min)
{
i -> TS = TIME;
i -> DIST = d_min + 1;
}
else
{
/* no parent is found */
add_to_changed_list(i);
/* process neighbors */
for (a0=i->first; a0; a0=a0->next)
{
j = a0 -> head;
if (j->is_sink && (a=j->parent))
{
if (a0->r_cap) set_active(j);
if (a!=TERMINAL && a!=ORPHAN && a->head==i)
{
set_orphan_rear(j); // add j to the end of the adoption list
}
}
}
}
}
/***********************************************************************/
template <typename captype, typename tcaptype, typename flowtype>
flowtype Graph<captype,tcaptype,flowtype>::maxflow(bool reuse_trees, Block<node_id>* _changed_list)
{
node *i, *j, *current_node = NULL;
arc *a;
nodeptr *np, *np_next;
if (!nodeptr_block)
{
nodeptr_block = new DBlock<nodeptr>(NODEPTR_BLOCK_SIZE, error_function);
}
changed_list = _changed_list;
if (maxflow_iteration == 0 && reuse_trees) { if (error_function) (*error_function)("reuse_trees cannot be used in the first call to maxflow()!"); exit(1); }
if (changed_list && !reuse_trees) { if (error_function) (*error_function)("changed_list cannot be used without reuse_trees!"); exit(1); }
if (reuse_trees) maxflow_reuse_trees_init();
else maxflow_init();
// main loop
while ( 1 )
{
// test_consistency(current_node);
if ((i=current_node))
{
i -> next = NULL; /* remove active flag */
if (!i->parent) i = NULL;
}
if (!i)
{
if (!(i = next_active())) break;
}
/* growth */
if (!i->is_sink)
{
/* grow source tree */
for (a=i->first; a; a=a->next)
if (a->r_cap)
{
j = a -> head;
if (!j->parent)
{
j -> is_sink = 0;
j -> parent = a -> sister;
j -> TS = i -> TS;
j -> DIST = i -> DIST + 1;
set_active(j);
add_to_changed_list(j);
}
else if (j->is_sink) break;
else if (j->TS <= i->TS &&
j->DIST > i->DIST)
{
/* heuristic - trying to make the distance from j to the source shorter */
j -> parent = a -> sister;
j -> TS = i -> TS;
j -> DIST = i -> DIST + 1;
}
}
}
else
{
/* grow sink tree */
for (a=i->first; a; a=a->next)
if (a->sister->r_cap)
{
j = a -> head;
if (!j->parent)
{
j -> is_sink = 1;
j -> parent = a -> sister;
j -> TS = i -> TS;
j -> DIST = i -> DIST + 1;
set_active(j);
add_to_changed_list(j);
}
else if (!j->is_sink) { a = a -> sister; break; }
else if (j->TS <= i->TS &&
j->DIST > i->DIST)
{
/* heuristic - trying to make the distance from j to the sink shorter */
j -> parent = a -> sister;
j -> TS = i -> TS;
j -> DIST = i -> DIST + 1;
}
}
}
TIME ++;
if (a)
{
i -> next = i; /* set active flag */
current_node = i;
/* augmentation */
augment(a);
/* augmentation end */
/* adoption */
while ((np=orphan_first))
{
np_next = np -> next;
np -> next = NULL;
while ((np=orphan_first))
{
orphan_first = np -> next;
i = np -> ptr;
nodeptr_block -> Delete(np);
if (!orphan_first) orphan_last = NULL;
if (i->is_sink) process_sink_orphan(i);
else process_source_orphan(i);
}
orphan_first = np_next;
}
/* adoption end */
}
else current_node = NULL;
}
// test_consistency();
if (!reuse_trees || (maxflow_iteration % 64) == 0)
{
delete nodeptr_block;
nodeptr_block = NULL;
}
maxflow_iteration ++;
return flow;
}
/***********************************************************************/
template <typename captype, typename tcaptype, typename flowtype>
void Graph<captype,tcaptype,flowtype>::test_consistency(node* current_node)
{
node *i;
arc *a;
int r;
int num1 = 0, num2 = 0;
// test whether all nodes i with i->next!=NULL are indeed in the queue
for (i=nodes; i<node_last; i++)
{
if (i->next || i==current_node) num1 ++;
}
for (r=0; r<3; r++)
{
i = (r == 2) ? current_node : queue_first[r];
if (i)
for ( ; ; i=i->next)
{
num2 ++;
if (i->next == i)
{
if (r<2) assert(i == queue_last[r]);
else assert(i == current_node);
break;
}
}
}
assert(num1 == num2);
for (i=nodes; i<node_last; i++)
{
// test whether all edges in seach trees are non-saturated
if (i->parent == NULL) {}
else if (i->parent == ORPHAN) {}
else if (i->parent == TERMINAL)
{
if (!i->is_sink) assert(i->tr_cap > 0);
else assert(i->tr_cap < 0);
}
else
{
if (!i->is_sink) assert (i->parent->sister->r_cap > 0);
else assert (i->parent->r_cap > 0);
}
// test whether passive nodes in search trees have neighbors in
// a different tree through non-saturated edges
if (i->parent && !i->next)
{
if (!i->is_sink)
{
assert(i->tr_cap >= 0);
for (a=i->first; a; a=a->next)
{
if (a->r_cap > 0) assert(a->head->parent && !a->head->is_sink);
}
}
else
{
assert(i->tr_cap <= 0);
for (a=i->first; a; a=a->next)
{
if (a->sister->r_cap > 0) assert(a->head->parent && a->head->is_sink);
}
}
}
// test marking invariants
if (i->parent && i->parent!=ORPHAN && i->parent!=TERMINAL)
{
assert(i->TS <= i->parent->head->TS);
if (i->TS == i->parent->head->TS) assert(i->DIST > i->parent->head->DIST);
}
}
}
#include "instances.inc"
/* block.h */
/*
Template classes Block and DBlock
Implement adding and deleting items of the same type in blocks.
If there there are many items then using Block or DBlock
is more efficient than using 'new' and 'delete' both in terms
of memory and time since
(1) On some systems there is some minimum amount of memory
that 'new' can allocate (e.g., 64), so if items are
small that a lot of memory is wasted.
(2) 'new' and 'delete' are designed for items of varying size.
If all items has the same size, then an algorithm for
adding and deleting can be made more efficient.
(3) All Block and DBlock functions are inline, so there are
no extra function calls.
Differences between Block and DBlock:
(1) DBlock allows both adding and deleting items,
whereas Block allows only adding items.
(2) Block has an additional operation of scanning
items added so far (in the order in which they were added).
(3) Block allows to allocate several consecutive
items at a time, whereas DBlock can add only a single item.
Note that no constructors or destructors are called for items.
Example usage for items of type 'MyType':
///
#include "block.h"
#define BLOCK_SIZE 1024
typedef struct { int a, b; } MyType;
MyType *ptr, *array[10000];
...
Block<MyType> *block = new Block<MyType>(BLOCK_SIZE);
// adding items
for (int i=0; i<sizeof(array); i++)
{
ptr = block -> New();
ptr -> a = ptr -> b = rand();
}
// reading items
for (ptr=block->ScanFirst(); ptr; ptr=block->ScanNext())
{
printf("%d %d\n", ptr->a, ptr->b);
}
delete block;
...
DBlock<MyType> *dblock = new DBlock<MyType>(BLOCK_SIZE);
// adding items
for (int i=0; i<sizeof(array); i++)
{
array[i] = dblock -> New();
}
// deleting items
for (int i=0; i<sizeof(array); i+=2)
{
dblock -> Delete(array[i]);
}
// adding items
for (int i=0; i<sizeof(array); i++)
{
array[i] = dblock -> New();
}
delete dblock;
///
Note that DBlock deletes items by marking them as
empty (i.e., by adding them to the list of free items),
so that this memory could be used for subsequently
added items. Thus, at each moment the memory allocated
is determined by the maximum number of items allocated
simultaneously at earlier moments. All memory is
deallocated only when the destructor is called.
*/
#ifndef __BLOCK_H__
#define __BLOCK_H__
#include <stdlib.h>
/***********************************************************************/
/***********************************************************************/
/***********************************************************************/
template <class Type> class Block
{
public:
/* Constructor. Arguments are the block size and
(optionally) the pointer to the function which
will be called if allocation failed; the message
passed to this function is "Not enough memory!" */
Block(int size, void (*err_function)(char *) = NULL) { first = last = NULL; block_size = size; error_function = err_function; }
/* Destructor. Deallocates all items added so far */
~Block() { while (first) { block *next = first -> next; delete[] ((char*)first); first = next; } }
/* Allocates 'num' consecutive items; returns pointer
to the first item. 'num' cannot be greater than the
block size since items must fit in one block */
Type *New(int num = 1)
{
Type *t;
if (!last || last->current + num > last->last)
{
if (last && last->next) last = last -> next;
else
{
block *next = (block *) new char [sizeof(block) + (block_size-1)*sizeof(Type)];
if (!next) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
if (last) last -> next = next;
else first = next;
last = next;
last -> current = & ( last -> data[0] );
last -> last = last -> current + block_size;
last -> next = NULL;
}
}
t = last -> current;
last -> current += num;
return t;
}
/* Returns the first item (or NULL, if no items were added) */
Type *ScanFirst()
{
for (scan_current_block=first; scan_current_block; scan_current_block = scan_current_block->next)
{
scan_current_data = & ( scan_current_block -> data[0] );
if (scan_current_data < scan_current_block -> current) return scan_current_data ++;
}
return NULL;
}
/* Returns the next item (or NULL, if all items have been read)
Can be called only if previous ScanFirst() or ScanNext()
call returned not NULL. */
Type *ScanNext()
{
while (scan_current_data >= scan_current_block -> current)
{
scan_current_block = scan_current_block -> next;
if (!scan_current_block) return NULL;
scan_current_data = & ( scan_current_block -> data[0] );
}
return scan_current_data ++;
}
/* Marks all elements as empty */
void Reset()
{
block *b;
if (!first) return;
for (b=first; ; b=b->next)
{
b -> current = & ( b -> data[0] );
if (b == last) break;
}
last = first;
}
/***********************************************************************/
private:
typedef struct block_st
{
Type *current, *last;
struct block_st *next;
Type data[1];
} block;
int block_size;
block *first;
block *last;
block *scan_current_block;
Type *scan_current_data;
void (*error_function)(char *);
};
/***********************************************************************/
/***********************************************************************/
/***********************************************************************/
template <class Type> class DBlock
{
public:
/* Constructor. Arguments are the block size and
(optionally) the pointer to the function which
will be called if allocation failed; the message
passed to this function is "Not enough memory!" */
DBlock(int size, void (*err_function)(char *) = NULL) { first = NULL; first_free = NULL; block_size = size; error_function = err_function; }
/* Destructor. Deallocates all items added so far */
~DBlock() { while (first) { block *next = first -> next; delete[] ((char*)first); first = next; } }
/* Allocates one item */
Type *New()
{
block_item *item;
if (!first_free)
{
block *next = first;
first = (block *) new char [sizeof(block) + (block_size-1)*sizeof(block_item)];
if (!first) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
first_free = & (first -> data[0] );
for (item=first_free; item<first_free+block_size-1; item++)
item -> next_free = item + 1;
item -> next_free = NULL;
first -> next = next;
}
item = first_free;
first_free = item -> next_free;
return (Type *) item;
}
/* Deletes an item allocated previously */
void Delete(Type *t)
{
((block_item *) t) -> next_free = first_free;
first_free = (block_item *) t;
}
/***********************************************************************/
private:
typedef union block_item_st
{
Type t;
block_item_st *next_free;
} block_item;
typedef struct block_st
{
struct block_st *next;
block_item data[1];
} block;
int block_size;
block *first;
block_item *first_free;
void (*error_function)(char *);
};
#endif
/* graph.h */
/*
This software library implements the maxflow algorithm
described in
"An Experimental Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision."
Yuri Boykov and Vladimir Kolmogorov.
In IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI),
September 2004
This algorithm was developed by Yuri Boykov and Vladimir Kolmogorov
at Siemens Corporate Research. To make it available for public use,
it was later reimplemented by Vladimir Kolmogorov based on open publications.
If you use this software for research purposes, you should cite
the aforementioned paper in any resulting publication.
----------------------------------------------------------------------
REUSING TREES:
Starting with version 3.0, there is a also an option of reusing search
trees from one maxflow computation to the next, as described in
"Efficiently Solving Dynamic Markov Random Fields Using Graph Cuts."
Pushmeet Kohli and Philip H.S. Torr
International Conference on Computer Vision (ICCV), 2005
If you use this option, you should cite
the aforementioned paper in any resulting publication.
*/
/*
For description, license, example usage see README.TXT.
*/
#ifndef __GRAPH_H__
#define __GRAPH_H__
#include <string.h>
#include "block.h"
#include <assert.h>
// NOTE: in UNIX you need to use -DNDEBUG preprocessor option to supress assert's!!!
// captype: type of edge capacities (excluding t-links)
// tcaptype: type of t-links (edges between nodes and terminals)
// flowtype: type of total flow
//
// Current instantiations are in instances.inc
template <typename captype, typename tcaptype, typename flowtype> class Graph
{
public:
typedef enum
{
SOURCE = 0,
SINK = 1
} termtype; // terminals
typedef int node_id;
/
// BASIC INTERFACE FUNCTIONS //
// (should be enough for most applications) //
/
// Constructor.
// The first argument gives an estimate of the maximum number of nodes that can be added
// to the graph, and the second argument is an estimate of the maximum number of edges.
// The last (optional) argument is the pointer to the function which will be called
// if an error occurs; an error message is passed to this function.
// If this argument is omitted, exit(1) will be called.
//
// IMPORTANT: It is possible to add more nodes to the graph than node_num_max
// (and node_num_max can be zero). However, if the count is exceeded, then
// the internal memory is reallocated (increased by 50%) which is expensive.
// Also, temporarily the amount of allocated memory would be more than twice than needed.
// Similarly for edges.
// If you wish to avoid this overhead, you can download version 2.2, where nodes and edges are stored in blocks.
Graph(int node_num_max, int edge_num_max, void (*err_function)(char *) = NULL);
// Destructor
~Graph();
// Adds node(s) to the graph. By default, one node is added (num=1); then first call returns 0, second call returns 1, and so on.
// If num>1, then several nodes are added, and node_id of the first one is returned.
// IMPORTANT: see note about the constructor
node_id add_node(int num = 1);
// Adds a bidirectional edge between 'i' and 'j' with the weights 'cap' and 'rev_cap'.
// IMPORTANT: see note about the constructor
void add_edge(node_id i, node_id j, captype cap, captype rev_cap);
// Adds new edges 'SOURCE->i' and 'i->SINK' with corresponding weights.
// Can be called multiple times for each node.
// Weights can be negative.
// NOTE: the number of such edges is not counted in edge_num_max.
// No internal memory is allocated by this call.
void add_tweights(node_id i, tcaptype cap_source, tcaptype cap_sink);
// Computes the maxflow. Can be called several times.
// FOR DESCRIPTION OF reuse_trees, SEE mark_node().
// FOR DESCRIPTION OF changed_list, SEE remove_from_changed_list().
flowtype maxflow(bool reuse_trees = false, Block<node_id>* changed_list = NULL);
// After the maxflow is computed, this function returns to which
// segment the node 'i' belongs (Graph<captype,tcaptype,flowtype>::SOURCE or Graph<captype,tcaptype,flowtype>::SINK).
//
// Occasionally there may be several minimum cuts. If a node can be assigned
// to both the source and the sink, then default_segm is returned.
termtype what_segment(node_id i, termtype default_segm = SOURCE);
//
// ADVANCED INTERFACE FUNCTIONS //
// (provide access to the graph) //
//
private:
struct node;
struct arc;
public:
// 1. Reallocating graph. //
// Removes all nodes and edges.
// After that functions add_node() and add_edge() must be called again.
//
// Advantage compared to deleting Graph and allocating it again:
// no calls to delete/new (which could be quite slow).
//
// If the graph structure stays the same, then an alternative
// is to go through all nodes/edges and set new residual capacities
// (see functions below).
void reset();
// 2. Functions for getting pointers to arcs and for reading graph structure. //
// NOTE: adding new arcs may invalidate these pointers (if reallocation //
// happens). So it's best not to add arcs while reading graph structure. //
// The following two functions return arcs in the same order that they
// were added to the graph. NOTE: for each call add_edge(i,j,cap,cap_rev)
// the first arc returned will be i->j, and the second j->i.
// If there are no more arcs, then the function can still be called, but
// the returned arc_id is undetermined.
typedef arc* arc_id;
arc_id get_first_arc();
arc_id get_next_arc(arc_id a);
// other functions for reading graph structure
int get_node_num() { return node_num; }
int get_arc_num() { return (int)(arc_last - arcs); }
void get_arc_ends(arc_id a, node_id& i, node_id& j); // returns i,j to that a = i->j
///
// 3. Functions for reading residual capacities. //
///
// returns residual capacity of SOURCE->i minus residual capacity of i->SINK
tcaptype get_trcap(node_id i);
// returns residual capacity of arc a
captype get_rcap(arc* a);
/
// 4. Functions for setting residual capacities. //
// NOTE: If these functions are used, the value of the flow //
// returned by maxflow() will not be valid! //
/
void set_trcap(node_id i, tcaptype trcap);
void set_rcap(arc* a, captype rcap);
// 5. Functions related to reusing trees & list of changed nodes. //
// If flag reuse_trees is true while calling maxflow(), then search trees
// are reused from previous maxflow computation.
// In this case before calling maxflow() the user must
// specify which parts of the graph have changed by calling mark_node():
// add_tweights(i),set_trcap(i) => call mark_node(i)
// add_edge(i,j),set_rcap(a) => call mark_node(i); mark_node(j)
//
// This option makes sense only if a small part of the graph is changed.
// The initialization procedure goes only through marked nodes then.
//
// mark_node(i) can either be called before or after graph modification.
// Can be called more than once per node, but calls after the first one
// do not have any effect.
//
// NOTE:
// - This option cannot be used in the first call to maxflow().
// - It is not necessary to call mark_node() if the change is ``not essential'',
// i.e. sign(trcap) is preserved for a node and zero/nonzero status is preserved for an arc.
// - To check that you marked all necessary nodes, you can call maxflow(false) after calling maxflow(true).
// If everything is correct, the two calls must return the same value of flow. (Useful for debugging).
void mark_node(node_id i);
// If changed_list is not NULL while calling maxflow(), then the algorithm
// keeps a list of nodes which could potentially have changed their segmentation label.
// Nodes which are not in the list are guaranteed to keep their old segmentation label (SOURCE or SINK).
// Example usage:
//
// typedef Graph<int,int,int> G;
// G* g = new Graph(nodeNum, edgeNum);
// Block<G::node_id>* changed_list = new Block<G::node_id>(128);
//
// ... // add nodes and edges
//
// g->maxflow(); // first call should be without arguments
// for (int iter=0; iter<10; iter++)
// {
// ... // change graph, call mark_node() accordingly
//
// g->maxflow(true, changed_list);
// G::node_id* ptr;
// for (ptr=changed_list->ScanFirst(); ptr; ptr=changed_list->ScanNext())
// {
// G::node_id i = *ptr; assert(i>=0 && i<nodeNum);
// g->remove_from_changed_list(i);
// // do something with node i...
// if (g->what_segment(i) == G::SOURCE) { ... }
// }
// changed_list->Reset();
// }
// delete changed_list;
//
// NOTE:
// - If changed_list option is used, then reuse_trees must be used as well.
// - In the example above, the user may omit calls g->remove_from_changed_list(i) and changed_list->Reset() in a given iteration.
// Then during the next call to maxflow(true, &changed_list) new nodes will be added to changed_list.
// - If the next call to maxflow() does not use option reuse_trees, then calling remove_from_changed_list()
// is not necessary. ("changed_list->Reset()" or "delete changed_list" should still be called, though).
void remove_from_changed_list(node_id i)
{
assert(i>=0 && i<node_num && nodes[i].is_in_changed_list);
nodes[i].is_in_changed_list = 0;
}
/
/
/
private:
// internal variables and functions
struct node
{
arc *first; // first outcoming arc
arc *parent; // node's parent
node *next; // pointer to the next active node
// (or to itself if it is the last node in the list)
int TS; // timestamp showing when DIST was computed
int DIST; // distance to the terminal
int is_sink : 1; // flag showing whether the node is in the source or in the sink tree (if parent!=NULL)
int is_marked : 1; // set by mark_node()
int is_in_changed_list : 1; // set by maxflow if
tcaptype tr_cap; // if tr_cap > 0 then tr_cap is residual capacity of the arc SOURCE->node
// otherwise -tr_cap is residual capacity of the arc node->SINK
};
struct arc
{
node *head; // node the arc points to
arc *next; // next arc with the same originating node
arc *sister; // reverse arc
captype r_cap; // residual capacity
};
struct nodeptr
{
node *ptr;
nodeptr *next;
};
static const int NODEPTR_BLOCK_SIZE = 128;
node *nodes, *node_last, *node_max; // node_last = nodes+node_num, node_max = nodes+node_num_max;
arc *arcs, *arc_last, *arc_max; // arc_last = arcs+2*edge_num, arc_max = arcs+2*edge_num_max;
int node_num;
DBlock<nodeptr> *nodeptr_block;
void (*error_function)(char *); // this function is called if a error occurs,
// with a corresponding error message
// (or exit(1) is called if it's NULL)
flowtype flow; // total flow
// reusing trees & list of changed pixels
int maxflow_iteration; // counter
Block<node_id> *changed_list;
/
node *queue_first[2], *queue_last[2]; // list of active nodes
nodeptr *orphan_first, *orphan_last; // list of pointers to orphans
int TIME; // monotonically increasing global counter
/
void reallocate_nodes(int num); // num is the number of new nodes
void reallocate_arcs();
// functions for processing active list
void set_active(node *i);
node *next_active();
// functions for processing orphans list
void set_orphan_front(node* i); // add to the beginning of the list
void set_orphan_rear(node* i); // add to the end of the list
void add_to_changed_list(node* i);
void maxflow_init(); // called if reuse_trees == false
void maxflow_reuse_trees_init(); // called if reuse_trees == true
void augment(arc *middle_arc);
void process_source_orphan(node *i);
void process_sink_orphan(node *i);
void test_consistency(node* current_node=NULL); // debug function
};
///
// Implementation - inline functions //
///
template <typename captype, typename tcaptype, typename flowtype>
inline typename Graph<captype,tcaptype,flowtype>::node_id Graph<captype,tcaptype,flowtype>::add_node(int num)
{
assert(num > 0);
if (node_last + num > node_max) reallocate_nodes(num);
if (num == 1)
{
node_last -> first = NULL;
node_last -> tr_cap = 0;
node_last -> is_marked = 0;
node_last -> is_in_changed_list = 0;
node_last ++;
return node_num ++;
}
else
{
memset(node_last, 0, num*sizeof(node));
node_id i = node_num;
node_num += num;
node_last += num;
return i;
}
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::add_tweights(node_id i, tcaptype cap_source, tcaptype cap_sink)
{
assert(i >= 0 && i < node_num);
tcaptype delta = nodes[i].tr_cap;
if (delta > 0) cap_source += delta;
else cap_sink -= delta;
flow += (cap_source < cap_sink) ? cap_source : cap_sink;
nodes[i].tr_cap = cap_source - cap_sink;
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::add_edge(node_id _i, node_id _j, captype cap, captype rev_cap)
{
assert(_i >= 0 && _i < node_num);
assert(_j >= 0 && _j < node_num);
assert(_i != _j);
assert(cap >= 0);
assert(rev_cap >= 0);
if (arc_last == arc_max) reallocate_arcs();
arc *a = arc_last ++;
arc *a_rev = arc_last ++;
node* i = nodes + _i;
node* j = nodes + _j;
a -> sister = a_rev;
a_rev -> sister = a;
a -> next = i -> first;
i -> first = a;
a_rev -> next = j -> first;
j -> first = a_rev;
a -> head = j;
a_rev -> head = i;
a -> r_cap = cap;
a_rev -> r_cap = rev_cap;
}
template <typename captype, typename tcaptype, typename flowtype>
inline typename Graph<captype,tcaptype,flowtype>::arc* Graph<captype,tcaptype,flowtype>::get_first_arc()
{
return arcs;
}
template <typename captype, typename tcaptype, typename flowtype>
inline typename Graph<captype,tcaptype,flowtype>::arc* Graph<captype,tcaptype,flowtype>::get_next_arc(arc* a)
{
return a + 1;
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::get_arc_ends(arc* a, node_id& i, node_id& j)
{
assert(a >= arcs && a < arc_last);
i = (node_id) (a->sister->head - nodes);
j = (node_id) (a->head - nodes);
}
template <typename captype, typename tcaptype, typename flowtype>
inline tcaptype Graph<captype,tcaptype,flowtype>::get_trcap(node_id i)
{
assert(i>=0 && i<node_num);
return nodes[i].tr_cap;
}
template <typename captype, typename tcaptype, typename flowtype>
inline captype Graph<captype,tcaptype,flowtype>::get_rcap(arc* a)
{
assert(a >= arcs && a < arc_last);
return a->r_cap;
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::set_trcap(node_id i, tcaptype trcap)
{
assert(i>=0 && i<node_num);
nodes[i].tr_cap = trcap;
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::set_rcap(arc* a, captype rcap)
{
assert(a >= arcs && a < arc_last);
a->r_cap = rcap;
}
template <typename captype, typename tcaptype, typename flowtype>
inline typename Graph<captype,tcaptype,flowtype>::termtype Graph<captype,tcaptype,flowtype>::what_segment(node_id i, termtype default_segm)
{
if (nodes[i].parent)
{
return (nodes[i].is_sink) ? SINK : SOURCE;
}
else
{
return default_segm;
}
}
template <typename captype, typename tcaptype, typename flowtype>
inline void Graph<captype,tcaptype,flowtype>::mark_node(node_id _i)
{
node* i = nodes + _i;
if (!i->next)
{
/* it's not in the list yet */
if (queue_last[1]) queue_last[1] -> next = i;
else queue_first[1] = i;
queue_last[1] = i;
i -> next = i;
}
i->is_marked = 1;
}
#endif
#include "graph.h"
#ifdef _MSC_VER
#pragma warning(disable: 4661)
#endif
// Instantiations: <captype, tcaptype, flowtype>
// IMPORTANT:
// flowtype should be 'larger' than tcaptype
// tcaptype should be 'larger' than captype
template class Graph<int,int,int>;
template class Graph<short,int,int>;
template class Graph<float,float,float>;
template class Graph<double,double,double>;
使用了maxflow库：http://vision.csd.uwo.ca/code/ 中Max-flow/min-cut