Grab Cut学习笔记1(new min-cut /max-flow algorithm)

时间:2022-02-28 06:43:09

参(chao)考(xi)的资料还是很多的,好好学习,天天向上嘛~

1.swsamleo大神:http://blog.csdn.net/swsamleo/article/details/7915316

2.wstcegg大神:http://blog.csdn.net/wstcegg/article/details/39495535

3.zouxy09大神图像分割系列:http://blog.csdn.net/zouxy09/article/details/8532106

4.GraphCuts on the GPU:http://blog.sina.com.cn/s/blog_60a0e97e0101bc75.html

5.GraphCut,Max-flow/min-cut等的代码:

http://vision.csd.uwo.ca/code/

http://pub.ist.ac.at/~vnk/software.html

论文:

1.An Experimental Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision -- Yuri Y.Boykov

作者写了一个new min-cut /max-flow algorithm的实现:http://vision.csd.uwo.ca/code/maxflow-v3.01.zip

2.Interactive Graph Cuts for Optimal Boundary & Region Segmentation of Objects in N-D Images -- Yuri Y. Boykov

3. grab cut-Interactive Foreground Extraction using Iterated Graph Cuts

下面是maxflow的讲解,这个maxflow算法的作者称做new min-cut /max-flow algorithm,他与标准maxflow算法不太一样,是通过建立树的形式去寻找s->t路径的,算法复杂度陡升啊...但是据说对于图像而言,这个方案反而比标准算法更优...

Example of the search trees S (red nodes) and T (blue nodes) at the end of the growth stage when a path (yellow line) from the source s to the sink t is found. 

Active and passive nodes are labeled by letters A and P, correspondingly. Free nodes appear in black.


 Grab Cut学习笔记1(new min-cut /max-flow algorithm)


S树和T树最前沿的点称为active node,这些点的任务是去发展新的node。而被active node包围起来的那些点,则称之为passive node。而没有被发掘出来的点则称之为free node

An augmenting path is found as soon as an active node in one of the trees detects a neighboring node that belongs to the other tree.

 

The algorithm iteratively repeats the following three stages:

1“growth” stage: search trees S and T grow until they touch giving an s ->t path, grow S & T search trees, find an edge connecting them

首先算法分别从sourcesink出发,建立两棵广度优先搜索树。In tree S, all edges from each parent node to its children are nonsaturated, while, in tree T , 

edges from children to their parents are nonsaturated,直到找到一条s->t的路径。

2 “augmentation” stage: the found path is augmented, search tree(s) break into forest(s)

The augmentation stage augments the path found at the growth stage. Since we push through the largest flow possible, 

some edge(s) in the path become saturated.

增长第一阶段找到的s->t的路径,因为我们会在这条路径中找最大的流,所以S树和T树上的一些节点会变成暂时的孤儿节点,

我们将这些节点加入孤立节点集合。 这一步会把S树和T树打乱,分成好多棵树。that is, the edges linking them to their parents are no longer valid (they are saturated). In fact, the augmentation phase may 

split the search trees S and T into forests. The source s and the sink t are still roots of two of the trees, 

while orphans form roots of all other trees.

orphan集合中在增广时建立,每次更新一个边后,如果发现改变后的边的残留流量=0,则把边指向的那个点加入orphan set

3 “adoption” stage: trees S and T are restored.

The goal of the adoption stage is to restore the single-tree structure of sets S and T with roots in the source and the sink. At this stage, we try to find a new valid parent for each orphan. 

A new parent should belong to the same set, S or T , as the orphan. A parent should also be connected through a nonsaturated edge. 

If there is no qualifying parent, we remove the orphan from S or T and make it a free node. We also declare all its former children orphans. 

The stage terminates when no orphans are left and, thus, the search tree structures of S and T are restored. 

Since some orphan nodes in S and T may become free, the adoption stage results in contraction of these sets.

adoption阶段,反复从orphan set中取点,每取出一个孤儿,首先看看能不能找到一个新的parent,如果找不到则令其变成free node

并把他的child变成orphan,加入到orphan集合中。循环往复,直到orphan set = 空集。

嗯,上面的英文都是论文里的讲述,还是蛮清晰的。。

 

找到一个可行流后,要进行augmention,然后必然会出现饱和的边。

举个栗子:比如这样的一条路中:s->p1->p2->p3->p4->t,增广后p2->p3这条边饱和了。如果不做处理,再次找到一条路后,

那么这条路中就有可能包含一些已经饱和路径,那就没法进行增广了。所以,我们必须要去调整那些饱和的边,

使得在已经构建的路径中不存在饱和的边。

在本例中,p3称之为orphan(孤点),那么接下来就得做一个adopt orphan的操作,

意思就是说给每个p3这样的孤儿点找一个新的parent,养育完就没有orphan存在了。

方法如下:

检查p3neighbors,看看有没有一个neighbor q

(1)q->p3满足容量大于0

(2)q是已经被搜过的点,也就是说q已经在我们的了,in another word,当我们沿着汇点逆流而上搜索到q时,可以顺利地找到qfather

从而最终可以顺藤摸瓜摸到source上。

(3)通过q最终能到达source或者sink。因为有可能顺着q走着走着,最终走到一个free node去了。或者走到一个orphan去了,

而这个orphan最终也无人抚养而变成free node

tip1:不必每次要追溯到source/sink才罢休。

因为对于orphan的每一个邻居进行判断其是否originate from TERMINATE node时,都要逆流追溯至source / sink点。

那么其中的一些点可能要被追溯多次,那么这是一个重复的操作。如果一个点,已经被证明了,他是可以到达source / sink的,

那么下次当有点经过他的时候,他就可以直接告诉该点,OK,哥们歇歇吧,你是valid的,通过我可以追溯到source/sink

这就是在algorithm tunning里提到的mark的意思。

tip2:选择离source/sink最近的那个neighbor,作为orphan的新parent.

要实现这个优化,必须给每个节点附加一个属性:该节点到source/sink的最短距离。

并且在growth阶段,当一个active node q1 遇到另一个active node q2时,要比较一下是否把parent(q2)=q1后,

q2source/sink的距离更短,如果是,那么就调整下q2,使得它变成q1child

大神说:人往高处走,水往低处流,节点都往终点凑!

 

如果找不到这样的q,那么p3就不得不变成free node。然后再做以下两个调整:

(4)对于p3的邻居pk,如果pkp3的边(pk-->p3)的容量大于0, 则将pk设为active

(5) 对于p3的邻居ph,如果p3=parent(ph),那么把ph也设置为orphan,加入到orphan集合里。

这是因为当有一条路从sink走到ph的时候,会发现无路可走了,因为parent(ph)是一个free node!

所以我们必须对ph也做相同的处理,要么找一个新的parent,要么您老人家变成free node,先一边凉快会!

 

Opencv也有个这个算法的实现,opencvgrabcut算法就是调用的介个算法求的最小割。感谢wstcegggcgraph.hpp注解:

/*M///////////////////////////////////////////////////////////////////////////////////////
//
// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
// By downloading, copying, installing or using the software you agree to this license.
// If you do not agree to this license, do not download, install,
// copy or use the software.
//
//
// Intel License Agreement
// For Open Source Computer Vision Library
//
// Copyright (C) 2000, Intel Corporation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// * Redistribution's of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
//
// * Redistribution's in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// * The name of Intel Corporation may not be used to endorse or promote products
// derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/

#ifndef _CV_GCGRAPH_H_
#define _CV_GCGRAPH_H_


template <class TWeight> class GCGraph
{
public:
GCGraph();
GCGraph( unsigned int vtxCount, unsigned int edgeCount );
~GCGraph();
void create( unsigned int vtxCount, unsigned int edgeCount );//建立有vtxCount个顶点edgeCount条边的图哇
int addVtx();//添加顶点,返回这个顶点在图中的索引
void addEdges( int i, int j, TWeight w, TWeight revw );//i->j权值w,j->i权值revw
void addTermWeights( int i, TWeight sourceW, TWeight sinkW );//s->i,i->t的权值分别为sourceW,sinkW
TWeight maxFlow();//求最大流
bool inSourceSegment( int i );
private:
class Vtx
{
public:
Vtx *next; // initialized and used in maxFlow() only
int parent;//父节点->Vtx边(父边)
int first;//顶点的第一条边
int ts;//timestamp showing when DIST was computed
int dist;//distance to the terminal
TWeight weight;/*if weight > 0 then t is residual capacity of the arc SOURCE->node
otherwise -weight is residual capacity of the arc node->SINK */
uchar t;// flag showing whether the node is in the source or in the sink tree
};
class Edge
{
public:
int dst;//边的目标
int next;//同一个出发点的下一条边
TWeight weight;
};

std::vector<Vtx> vtcs;//顶点集合
std::vector<Edge> edges;//边的集合
TWeight flow;//图的流量
};

template <class TWeight>
GCGraph<TWeight>::GCGraph()
{
flow = 0;
}
template <class TWeight>
GCGraph<TWeight>::GCGraph( unsigned int vtxCount, unsigned int edgeCount )
{
create( vtxCount, edgeCount );
}
template <class TWeight>
GCGraph<TWeight>::~GCGraph()
{
}
template <class TWeight>
void GCGraph<TWeight>::create( unsigned int vtxCount, unsigned int edgeCount )
{
vtcs.reserve( vtxCount );
edges.reserve( edgeCount + 2 );//加两条到原点和汇点的边
flow = 0;
}

template <class TWeight>
int GCGraph<TWeight>::addVtx()
{
Vtx v;
memset( &v, 0, sizeof(Vtx));
vtcs.push_back(v);
return (int)vtcs.size() - 1;//返回顶点编号
}

template <class TWeight>
void GCGraph<TWeight>::addEdges( int i, int j, TWeight w, TWeight revw )
{
CV_Assert( i>=0 && i<(int)vtcs.size() );
CV_Assert( j>=0 && j<(int)vtcs.size() );
CV_Assert( w>=0 && revw>=0 );
CV_Assert( i != j );

if( !edges.size() )
edges.resize( 2 );//初始为至少两条边,即s到点,点到t的边

Edge fromI, toI;
fromI.dst = j;//边i->j的目标

//每个顶点出发的边(4个方向)建立一个链表,头插法插入
fromI.next = vtcs[i].first;//i出发的下一条边
fromI.weight = w;
vtcs[i].first = (int)edges.size();//顶点i的第一条边修改为当前的边
edges.push_back( fromI );

toI.dst = i;
toI.next = vtcs[j].first;
toI.weight = revw;
vtcs[j].first = (int)edges.size();
edges.push_back( toI );
}

template <class TWeight>
void GCGraph<TWeight>::addTermWeights( int i, TWeight sourceW, TWeight sinkW )
{
CV_Assert( i>=0 && i<(int)vtcs.size() );

TWeight dw = vtcs[i].weight;
if( dw > 0 )
sourceW += dw;
else
sinkW -= dw;
flow += (sourceW < sinkW) ? sourceW : sinkW;//流的大小
vtcs[i].weight = sourceW - sinkW;//顶点i的t-vaule初始为sourceW - sinkW;
}

template <class TWeight>
TWeight GCGraph<TWeight>::maxFlow()
{
const int TERMINAL = -1/*终节点*/, ORPHAN = -2/*孤立点*/;
Vtx stub, *nilNode = &stub, *first = nilNode, *last = nilNode;//队列保存当前活动结点,stub为哨兵结点
int curr_ts = 0;//当前时间
stub.next = nilNode;//初始化队列
Vtx *vtxPtr = &vtcs[0];
Edge *edgePtr = &edges[0];

std::vector<Vtx*> orphans;//孤立点集合,例如v->u边的权值为0,则把u加入孤立点集合

// initialize the active queue and the graph vertices
// 初始化活动结点(active node)队列,这些点的任务是去发展新的node
for( int i = 0; i < (int)vtcs.size(); i++ )
{
Vtx* v = vtxPtr + i;
v->ts = 0;
if( v->weight != 0 )
{
last = last->next = v;
v->dist = 1;
v->parent = TERMINAL;//初始其父边为-1
v->t = v->weight < 0;//v->weight < 0,则v为背景,否则为前景
}
else
v->parent = 0;//初始其父边为0,定义为孤立点,其到原点和汇点的权值相同
}
first = first->next;
last->next = nilNode;
nilNode->next = 0;

// run the search-path -> augment-graph -> restore-trees loop
// BFS搜索路径->增广->树的重构
for(;;)
{
Vtx* v, *u;//u为v的相邻顶点
int e0 = -1, ei = 0, ej = 0;
TWeight minWeight, weight;//minWeight路径最小割(流量), weight当前流量
uchar vt;//前景背景点标记

// grow S & T search trees, find an edge connecting them
// growth S 和 T 树,找到一条s->t的路径
while( first != nilNode )//活动节点去去发展新的活动节点
{
v = first;
if( v->parent )
{
vt = v->t;
//广度优先搜索来搜索增广路径
for( ei = v->first; ei != 0; ei = edgePtr[ei].next )
{
//ei^vt:一个整数与0异或结果不变,与1异或,如果是奇数与1异或结果是奇数-1,偶数与1异或结果是偶数+1
//而我们这里建图的时候,toI = fromI+1,所以当前的边是toI,那么与1异或后就变成了对应的fromI,反之亦然

if( edgePtr[ei^vt].weight == 0 )//假如边的剩余流量为0,继续遍历请寻找汇点
continue;

u = vtxPtr+edgePtr[ei].dst;//v->u边的终点u

if( !u->parent )//adopt orphan
{
u->t = vt;
u->parent = ei ^ 1;//ei为toU,则边fromu为其父边
u->ts = v->ts;
u->dist = v->dist + 1;//说明s->t的路径经过了v->u
if( !u->next )
{
u->next = nilNode;
last = last->next = u;//u变为活动节点
}
continue;
}

if( u->t != vt )///u的标记与v的标记不同,表示找到一条s->t的路径辣
{
e0 = ei ^ vt;///e0记下s树与t树的节点相连接的辣条边
break;
}

//u的路径长度大于v的长度+1,且u的时间戳较早,即之前经过过u
if( u->dist > v->dist+1 && u->ts <= v->ts )
{
// reassign the parent
u->parent = ei ^ 1;
u->ts = v->ts;
u->dist = v->dist + 1;
}
}
if( e0 > 0 )///找到一条s->t的路径
break;
}
// exclude the vertex from the active list
first = first->next;
v->next = 0;
}

if( e0 <= 0 )///找不到s->t的路径,结束
break;

// find the minimum edge weight along the path
//查找路径中的最小权值

minWeight = edgePtr[e0].weight;
assert( minWeight > 0 );

// k = 1: source tree, k = 0: destination tree
// k=1: 回溯s树,k=0: 回溯t树
for( int k = 1; k >= 0; k-- )
{
for( v = vtxPtr+edgePtr[e0^k].dst;; v = vtxPtr+edgePtr[ei].dst )
{
if( (ei = v->parent) < 0 )//ei纪录当前点的父边,回溯至终端结点,退出
break;
weight = edgePtr[ei^k].weight;
minWeight = MIN(minWeight, weight);
assert( minWeight > 0 );
}
weight = fabs(v->weight);
minWeight = MIN(minWeight, weight);
assert( minWeight > 0 );
}

// modify weights of the edges along the path and collect orphans
edgePtr[e0].weight -= minWeight;
edgePtr[e0^1].weight += minWeight;
flow += minWeight;//当前流量

// k = 1: source tree, k = 0: destination tree
/*
augmentation,orphan集合中在增广时建立,每次更新一个边后,
如果发现改变后的边的残留流量=0则把边指向的那个点加入orphan set
*/
for( int k = 1; k >= 0; k-- )
{
for( v = vtxPtr+edgePtr[e0^k].dst;; v = vtxPtr+edgePtr[ei].dst )
{
if( (ei = v->parent) < 0 )
break;
edgePtr[ei^(k^1)].weight += minWeight;
if( (edgePtr[ei^k].weight -= minWeight) == 0 )
{
orphans.push_back(v);
v->parent = ORPHAN;
}
}

v->weight = v->weight + minWeight*(1-k*2);
if( v->weight == 0 )
{
orphans.push_back(v);
v->parent = ORPHAN;
}
}

// restore the search trees by finding new parents for the orphans
/*
adoption,反复从orphan set中取点,每取出一个孤儿,
看看能不能找到一个新的parent,选择离source/sink最近的那个neighbor,作为orphan的新parent.
如果找不到则令其变成free node,并把他的child变成orphan,加入到orphan集合中,重复直到orphan set = null
注意:orphans and free nodes have no parents。
*/
curr_ts++;
while( !orphans.empty() )
{
Vtx* v2 = orphans.back();
orphans.pop_back();

int d, minDist = INT_MAX;
e0 = 0;
vt = v2->t;

for( ei = v2->first; ei != 0; ei = edgePtr[ei].next )
{
if( edgePtr[ei^(vt^1)].weight == 0 )
continue;
u = vtxPtr+edgePtr[ei].dst;//v2->u

if( u->t != vt || u->parent == 0 )//顶点的标记不同或找到的顶点也是孤立点
continue;

// compute the distance to the tree root
for( d = 0;; )
{
if( u->ts == curr_ts )
{
d += u->dist;
break;
}
ej = u->parent;
d++;
if( ej < 0 )// TERMINAL = -1, ORPHAN = -2
{
if( ej == ORPHAN )//孤立点的距离无穷大
d = INT_MAX-1;
else
{
u->ts = curr_ts;
u->dist = 1;
}
break;
}
u = vtxPtr+edgePtr[ej].dst;
}

// update the distance
if( ++d < INT_MAX )
{
if( d < minDist )
{
minDist = d;
e0 = ei;
}
for( u = vtxPtr+edgePtr[ei].dst; u->ts != curr_ts; u = vtxPtr+edgePtr[u->parent].dst )
{
u->ts = curr_ts;
u->dist = --d;
}
}
}

if( (v2->parent = e0) > 0 )
{
v2->ts = curr_ts;
v2->dist = minDist;//更新最小距离
continue;
}

/* no parent is found ,处理其相邻顶点*/
v2->ts = 0;
for( ei = v2->first; ei != 0; ei = edgePtr[ei].next )
{
u = vtxPtr+edgePtr[ei].dst;//u为v2相邻顶点,ei为v2->u
ej = u->parent;
if( u->t != vt || !ej )//标记不一样
continue;

if( edgePtr[ei^(vt^1)].weight && !u->next )//如果u->v2的权值大于0, 则将u设为active
{
u->next = nilNode;
last = last->next = u;
}
if( ej > 0 && vtxPtr+edgePtr[ej].dst == v2 )//如果v2为u的父亲,则u设置为孤立点
{
orphans.push_back(u);
u->parent = ORPHAN;
}
}
}
}
return flow;
}

template <class TWeight>
bool GCGraph<TWeight>::inSourceSegment( int i )
{
CV_Assert( i>=0 && i<(int)vtcs.size() );
return vtcs[i].t == 0;//返回顶点i是否为前景
}

#endif