Grab Cut学习理解之（1）new min-cut /max-flow algorithm

/*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