gcgraph.hpp

00001 /*M///////////////////////////////////////////////////////////////////////////////////////
00002 //
00003 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
00004 //
00005 //  By downloading, copying, installing or using the software you agree to this license.
00006 //  If you do not agree to this license, do not download, install,
00007 //  copy or use the software.
00008 //
00009 //
00010 //                        Intel License Agreement
00011 //                For Open Source Computer Vision Library
00012 //
00013 // Copyright (C) 2000, Intel Corporation, all rights reserved.
00014 // Third party copyrights are property of their respective owners.
00015 //
00016 // Redistribution and use in source and binary forms, with or without modification,
00017 // are permitted provided that the following conditions are met:
00018 //
00019 //   * Redistribution's of source code must retain the above copyright notice,
00020 //     this list of conditions and the following disclaimer.
00021 //
00022 //   * Redistribution's in binary form must reproduce the above copyright notice,
00023 //     this list of conditions and the following disclaimer in the documentation
00024 //     and/or other materials provided with the distribution.
00025 //
00026 //   * The name of Intel Corporation may not be used to endorse or promote products
00027 //     derived from this software without specific prior written permission.
00028 //
00029 // This software is provided by the copyright holders and contributors "as is" and
00030 // any express or implied warranties, including, but not limited to, the implied
00031 // warranties of merchantability and fitness for a particular purpose are disclaimed.
00032 // In no event shall the Intel Corporation or contributors be liable for any direct,
00033 // indirect, incidental, special, exemplary, or consequential damages
00034 // (including, but not limited to, procurement of substitute goods or services;
00035 // loss of use, data, or profits; or business interruption) however caused
00036 // and on any theory of liability, whether in contract, strict liability,
00037 // or tort (including negligence or otherwise) arising in any way out of
00038 // the use of this software, even if advised of the possibility of such damage.
00039 //
00040 //M*/
00041 
00042 #ifndef _CV_GCGRAPH_H_
00043 #define _CV_GCGRAPH_H_
00044 
00045 template <class TWeight> class GCGraph
00046 {
00047 public:
00048     GCGraph();
00049     GCGraph( unsigned int vtxCount, unsigned int edgeCount );
00050     ~GCGraph();
00051     void create( unsigned int vtxCount, unsigned int edgeCount );
00052     int addVtx();
00053     void addEdges( int i, int j, TWeight w, TWeight revw );
00054     void addTermWeights( int i, TWeight sourceW, TWeight sinkW );
00055     TWeight maxFlow();
00056     bool inSourceSegment( int i );
00057 private:
00058     class Vtx
00059     {
00060     public:
00061         Vtx *next; // initialized and used in maxFlow() only
00062         int parent;
00063         int first;
00064         int ts;
00065         int dist;
00066         TWeight weight;
00067         uchar t;
00068     };
00069     class Edge
00070     {
00071     public:
00072         int dst;
00073         int next;
00074         TWeight weight;
00075     };
00076 
00077     std::vector<Vtx> vtcs;
00078     std::vector<Edge> edges;
00079     TWeight flow;
00080 };
00081 
00082 template <class TWeight>
00083 GCGraph<TWeight>::GCGraph()
00084 {
00085     flow = 0;
00086 }
00087 template <class TWeight>
00088 GCGraph<TWeight>::GCGraph( unsigned int vtxCount, unsigned int edgeCount )
00089 {
00090     create( vtxCount, edgeCount );
00091 }
00092 template <class TWeight>
00093 GCGraph<TWeight>::~GCGraph()
00094 {
00095 }
00096 template <class TWeight>
00097 void GCGraph<TWeight>::create( unsigned int vtxCount, unsigned int edgeCount )
00098 {
00099     vtcs.reserve( vtxCount );
00100     edges.reserve( edgeCount + 2 );
00101     flow = 0;
00102 }
00103 
00104 template <class TWeight>
00105 int GCGraph<TWeight>::addVtx()
00106 {
00107     Vtx v;
00108     memset( &v, 0, sizeof(Vtx));
00109     vtcs.push_back(v);
00110     return (int)vtcs.size() - 1;
00111 }
00112 
00113 template <class TWeight>
00114 void GCGraph<TWeight>::addEdges( int i, int j, TWeight w, TWeight revw )
00115 {
00116     CV_Assert( i>=0 && i<(int)vtcs.size() );
00117     CV_Assert( j>=0 && j<(int)vtcs.size() );
00118     CV_Assert( w>=0 && revw>=0 );
00119     CV_Assert( i != j );
00120 
00121     if( !edges.size() )
00122         edges.resize( 2 );
00123 
00124     Edge fromI, toI;
00125     fromI.dst = j;
00126     fromI.next = vtcs[i].first;
00127     fromI.weight = w;
00128     vtcs[i].first = (int)edges.size();
00129     edges.push_back( fromI );
00130 
00131     toI.dst = i;
00132     toI.next = vtcs[j].first;
00133     toI.weight = revw;
00134     vtcs[j].first = (int)edges.size();
00135     edges.push_back( toI );
00136 }
00137 
00138 template <class TWeight>
00139 void GCGraph<TWeight>::addTermWeights( int i, TWeight sourceW, TWeight sinkW )
00140 {
00141     CV_Assert( i>=0 && i<(int)vtcs.size() );
00142 
00143     TWeight dw = vtcs[i].weight;
00144     if( dw > 0 )
00145         sourceW += dw;
00146     else
00147         sinkW -= dw;
00148     flow += (sourceW < sinkW) ? sourceW : sinkW;
00149     vtcs[i].weight = sourceW - sinkW;
00150 }
00151 
00152 template <class TWeight>
00153 TWeight GCGraph<TWeight>::maxFlow()
00154 {
00155     const int TERMINAL = -1, ORPHAN = -2;
00156     Vtx stub, *nilNode = &stub, *first = nilNode, *last = nilNode;
00157     int curr_ts = 0;
00158     stub.next = nilNode;
00159     Vtx *vtxPtr = &vtcs[0];
00160     Edge *edgePtr = &edges[0];
00161 
00162     std::vector<Vtx*> orphans;
00163 
00164     // initialize the active queue and the graph vertices
00165     for( int i = 0; i < (int)vtcs.size(); i++ )
00166     {
00167         Vtx* v = vtxPtr + i;
00168         v->ts = 0;
00169         if( v->weight != 0 )
00170         {
00171             last = last->next = v;
00172             v->dist = 1;
00173             v->parent = TERMINAL;
00174             v->t = v->weight < 0;
00175         }
00176         else
00177             v->parent = 0;
00178     }
00179     first = first->next;
00180     last->next = nilNode;
00181     nilNode->next = 0;
00182 
00183     // run the search-path -> augment-graph -> restore-trees loop
00184     for(;;)
00185     {
00186         Vtx* v, *u;
00187         int e0 = -1, ei = 0, ej = 0;
00188         TWeight minWeight, weight;
00189         uchar vt;
00190 
00191         // grow S & T search trees, find an edge connecting them
00192         while( first != nilNode )
00193         {
00194             v = first;
00195             if( v->parent )
00196             {
00197                 vt = v->t;
00198                 for( ei = v->first; ei != 0; ei = edgePtr[ei].next )
00199                 {
00200                     if( edgePtr[ei^vt].weight == 0 )
00201                         continue;
00202                     u = vtxPtr+edgePtr[ei].dst;
00203                     if( !u->parent )
00204                     {
00205                         u->t = vt;
00206                         u->parent = ei ^ 1;
00207                         u->ts = v->ts;
00208                         u->dist = v->dist + 1;
00209                         if( !u->next )
00210                         {
00211                             u->next = nilNode;
00212                             last = last->next = u;
00213                         }
00214                         continue;
00215                     }
00216 
00217                     if( u->t != vt )
00218                     {
00219                         e0 = ei ^ vt;
00220                         break;
00221                     }
00222 
00223                     if( u->dist > v->dist+1 && u->ts <= v->ts )
00224                     {
00225                         // reassign the parent
00226                         u->parent = ei ^ 1;
00227                         u->ts = v->ts;
00228                         u->dist = v->dist + 1;
00229                     }
00230                 }
00231                 if( e0 > 0 )
00232                     break;
00233             }
00234             // exclude the vertex from the active list
00235             first = first->next;
00236             v->next = 0;
00237         }
00238 
00239         if( e0 <= 0 )
00240             break;
00241 
00242         // find the minimum edge weight along the path
00243         minWeight = edgePtr[e0].weight;
00244         CV_Assert( minWeight > 0 );
00245         // k = 1: source tree, k = 0: destination tree
00246         for( int k = 1; k >= 0; k-- )
00247         {
00248             for( v = vtxPtr+edgePtr[e0^k].dst;; v = vtxPtr+edgePtr[ei].dst )
00249             {
00250                 if( (ei = v->parent) < 0 )
00251                     break;
00252                 weight = edgePtr[ei^k].weight;
00253                 minWeight = MIN(minWeight, weight);
00254                 CV_Assert( minWeight > 0 );
00255             }
00256             weight = fabs(v->weight);
00257             minWeight = MIN(minWeight, weight);
00258             CV_Assert( minWeight > 0 );
00259         }
00260 
00261         // modify weights of the edges along the path and collect orphans
00262         edgePtr[e0].weight -= minWeight;
00263         edgePtr[e0^1].weight += minWeight;
00264         flow += minWeight;
00265 
00266         // k = 1: source tree, k = 0: destination tree
00267         for( int k = 1; k >= 0; k-- )
00268         {
00269             for( v = vtxPtr+edgePtr[e0^k].dst;; v = vtxPtr+edgePtr[ei].dst )
00270             {
00271                 if( (ei = v->parent) < 0 )
00272                     break;
00273                 edgePtr[ei^(k^1)].weight += minWeight;
00274                 if( (edgePtr[ei^k].weight -= minWeight) == 0 )
00275                 {
00276                     orphans.push_back(v);
00277                     v->parent = ORPHAN;
00278                 }
00279             }
00280 
00281             v->weight = v->weight + minWeight*(1-k*2);
00282             if( v->weight == 0 )
00283             {
00284                orphans.push_back(v);
00285                v->parent = ORPHAN;
00286             }
00287         }
00288 
00289         // restore the search trees by finding new parents for the orphans
00290         curr_ts++;
00291         while( !orphans.empty() )
00292         {
00293             Vtx* v2 = orphans.back();
00294             orphans.pop_back();
00295 
00296             int d, minDist = INT_MAX;
00297             e0 = 0;
00298             vt = v2->t;
00299 
00300             for( ei = v2->first; ei != 0; ei = edgePtr[ei].next )
00301             {
00302                 if( edgePtr[ei^(vt^1)].weight == 0 )
00303                     continue;
00304                 u = vtxPtr+edgePtr[ei].dst;
00305                 if( u->t != vt || u->parent == 0 )
00306                     continue;
00307                 // compute the distance to the tree root
00308                 for( d = 0;; )
00309                 {
00310                     if( u->ts == curr_ts )
00311                     {
00312                         d += u->dist;
00313                         break;
00314                     }
00315                     ej = u->parent;
00316                     d++;
00317                     if( ej < 0 )
00318                     {
00319                         if( ej == ORPHAN )
00320                             d = INT_MAX-1;
00321                         else
00322                         {
00323                             u->ts = curr_ts;
00324                             u->dist = 1;
00325                         }
00326                         break;
00327                     }
00328                     u = vtxPtr+edgePtr[ej].dst;
00329                 }
00330 
00331                 // update the distance
00332                 if( ++d < INT_MAX )
00333                 {
00334                     if( d < minDist )
00335                     {
00336                         minDist = d;
00337                         e0 = ei;
00338                     }
00339                     for( u = vtxPtr+edgePtr[ei].dst; u->ts != curr_ts; u = vtxPtr+edgePtr[u->parent].dst )
00340                     {
00341                         u->ts = curr_ts;
00342                         u->dist = --d;
00343                     }
00344                 }
00345             }
00346 
00347             if( (v2->parent = e0) > 0 )
00348             {
00349                 v2->ts = curr_ts;
00350                 v2->dist = minDist;
00351                 continue;
00352             }
00353 
00354             /* no parent is found */
00355             v2->ts = 0;
00356             for( ei = v2->first; ei != 0; ei = edgePtr[ei].next )
00357             {
00358                 u = vtxPtr+edgePtr[ei].dst;
00359                 ej = u->parent;
00360                 if( u->t != vt || !ej )
00361                     continue;
00362                 if( edgePtr[ei^(vt^1)].weight && !u->next )
00363                 {
00364                     u->next = nilNode;
00365                     last = last->next = u;
00366                 }
00367                 if( ej > 0 && vtxPtr+edgePtr[ej].dst == v2 )
00368                 {
00369                     orphans.push_back(u);
00370                     u->parent = ORPHAN;
00371                 }
00372             }
00373         }
00374     }
00375     return flow;
00376 }
00377 
00378 template <class TWeight>
00379 bool GCGraph<TWeight>::inSourceSegment( int i )
00380 {
00381     CV_Assert( i>=0 && i<(int)vtcs.size() );
00382     return vtcs[i].t == 0;
00383 }
00384 
00385 #endif
00386