modules/imgproc/src/gcgraph.hpp

   1 /*M///////////////////////////////////////////////////////////////////////////////////////
   2 //
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
   4 //
   5 //  By downloading, copying, installing or using the software you agree to this license.
   6 //  If you do not agree to this license, do not download, install,
   7 //  copy or use the software.
   8 //
   9 //
  10 //                        Intel License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
  14 // Third party copyrights are property of their respective owners.
  15 //
  16 // Redistribution and use in source and binary forms, with or without modification,
  17 // are permitted provided that the following conditions are met:
  18 //
  19 //   * Redistribution's of source code must retain the above copyright notice,
  20 //     this list of conditions and the following disclaimer.
  21 //
  22 //   * Redistribution's in binary form must reproduce the above copyright notice,
  23 //     this list of conditions and the following disclaimer in the documentation
  24 //     and/or other materials provided with the distribution.
  25 //
  26 //   * The name of Intel Corporation may not be used to endorse or promote products
  27 //     derived from this software without specific prior written permission.
  28 //
  29 // This software is provided by the copyright holders and contributors "as is" and
  30 // any express or implied warranties, including, but not limited to, the implied
  31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  32 // In no event shall the Intel Corporation or contributors be liable for any direct,
  33 // indirect, incidental, special, exemplary, or consequential damages
  34 // (including, but not limited to, procurement of substitute goods or services;
  35 // loss of use, data, or profits; or business interruption) however caused
  36 // and on any theory of liability, whether in contract, strict liability,
  37 // or tort (including negligence or otherwise) arising in any way out of
  38 // the use of this software, even if advised of the possibility of such damage.
  39 //
  40 //M*/
  41
  42 #ifndef _CV_GCGRAPH_H_
  43 #define _CV_GCGRAPH_H_
  44
  45 template <class TWeight> class GCGraph
  46 {
  47 public:
  48     GCGraph();
  49     GCGraph( unsigned int vtxCount, unsigned int edgeCount );
  50     ~GCGraph();
  51     void create( unsigned int vtxCount, unsigned int edgeCount );
  52     int addVtx();
  53     void addEdges( int i, int j, TWeight w, TWeight revw );
  54     void addTermWeights( int i, TWeight sourceW, TWeight sinkW );
  55     TWeight maxFlow();
  56     bool inSourceSegment( int i );
  57 private:
  58     class Vtx
  59     {
  60     public:
  61         Vtx *next; // initialized and used in maxFlow() only
  62         int parent;
  63         int first;
  64         int ts;
  65         int dist;
  66         TWeight weight;
  67         uchar t;
  68     };
  69     class Edge
  70     {
  71     public:
  72         int dst;
  73         int next;
  74         TWeight weight;
  75     };
  76
  77     std::vector<Vtx> vtcs;
  78     std::vector<Edge> edges;
  79     TWeight flow;
  80 };
  81
  82 template <class TWeight>
  83 GCGraph<TWeight>::GCGraph()
  84 {
  85     flow = 0;
  86 }
  87 template <class TWeight>
  88 GCGraph<TWeight>::GCGraph( unsigned int vtxCount, unsigned int edgeCount )
  89 {
  90     create( vtxCount, edgeCount );
  91 }
  92 template <class TWeight>
  93 GCGraph<TWeight>::~GCGraph()
  94 {
  95 }
  96 template <class TWeight>
  97 void GCGraph<TWeight>::create( unsigned int vtxCount, unsigned int edgeCount )
  98 {
  99     vtcs.reserve( vtxCount );
 100     edges.reserve( edgeCount + 2 );
 101     flow = 0;
 102 }
 103
 104 template <class TWeight>
 105 int GCGraph<TWeight>::addVtx()
 106 {
 107     Vtx v;
 108     memset( &v, 0, sizeof(Vtx));
 109     vtcs.push_back(v);
 110     return (int)vtcs.size() - 1;
 111 }
 112
 113 template <class TWeight>
 114 void GCGraph<TWeight>::addEdges( int i, int j, TWeight w, TWeight revw )
 115 {
 116     CV_Assert( i>=0 && i<(int)vtcs.size() );
 117     CV_Assert( j>=0 && j<(int)vtcs.size() );
 118     CV_Assert( w>=0 && revw>=0 );
 119     CV_Assert( i != j );
 120
 121     if( !edges.size() )
 122         edges.resize( 2 );
 123
 124     Edge fromI, toI;
 125     fromI.dst = j;
 126     fromI.next = vtcs[i].first;
 127     fromI.weight = w;
 128     vtcs[i].first = (int)edges.size();
 129     edges.push_back( fromI );
 130
 131     toI.dst = i;
 132     toI.next = vtcs[j].first;
 133     toI.weight = revw;
 134     vtcs[j].first = (int)edges.size();
 135     edges.push_back( toI );
 136 }
 137
 138 template <class TWeight>
 139 void GCGraph<TWeight>::addTermWeights( int i, TWeight sourceW, TWeight sinkW )
 140 {
 141     CV_Assert( i>=0 && i<(int)vtcs.size() );
 142
 143     TWeight dw = vtcs[i].weight;
 144     if( dw > 0 )
 145         sourceW += dw;
 146     else
 147         sinkW -= dw;
 148     flow += (sourceW < sinkW) ? sourceW : sinkW;
 149     vtcs[i].weight = sourceW - sinkW;
 150 }
 151
 152 template <class TWeight>
 153 TWeight GCGraph<TWeight>::maxFlow()
 154 {
 155     const int TERMINAL = -1, ORPHAN = -2;
 156     Vtx stub, *nilNode = &stub, *first = nilNode, *last = nilNode;
 157     int curr_ts = 0;
 158     stub.next = nilNode;
 159     Vtx *vtxPtr = &vtcs[0];
 160     Edge *edgePtr = &edges[0];
 161
 162     std::vector<Vtx*> orphans;
 163
 164     // initialize the active queue and the graph vertices
 165     for( int i = 0; i < (int)vtcs.size(); i++ )
 166     {
 167         Vtx* v = vtxPtr + i;
 168         v->ts = 0;
 169         if( v->weight != 0 )
 170         {
 171             last = last->next = v;
 172             v->dist = 1;
 173             v->parent = TERMINAL;
 174             v->t = v->weight < 0;
 175         }
 176         else
 177             v->parent = 0;
 178     }
 179     first = first->next;
 180     last->next = nilNode;
 181     nilNode->next = 0;
 182
 183     // run the search-path -> augment-graph -> restore-trees loop
 184     for(;;)
 185     {
 186         Vtx* v, *u;
 187         int e0 = -1, ei = 0, ej = 0;
 188         TWeight minWeight, weight;
 189         uchar vt;
 190
 191         // grow S & T search trees, find an edge connecting them
 192         while( first != nilNode )
 193         {
 194             v = first;
 195             if( v->parent )
 196             {
 197                 vt = v->t;
 198                 for( ei = v->first; ei != 0; ei = edgePtr[ei].next )
 199                 {
 200                     if( edgePtr[ei^vt].weight == 0 )
 201                         continue;
 202                     u = vtxPtr+edgePtr[ei].dst;
 203                     if( !u->parent )
 204                     {
 205                         u->t = vt;
 206                         u->parent = ei ^ 1;
 207                         u->ts = v->ts;
 208                         u->dist = v->dist + 1;
 209                         if( !u->next )
 210                         {
 211                             u->next = nilNode;
 212                             last = last->next = u;
 213                         }
 214                         continue;
 215                     }
 216
 217                     if( u->t != vt )
 218                     {
 219                         e0 = ei ^ vt;
 220                         break;
 221                     }
 222
 223                     if( u->dist > v->dist+1 && u->ts <= v->ts )
 224                     {
 225                         // reassign the parent
 226                         u->parent = ei ^ 1;
 227                         u->ts = v->ts;
 228                         u->dist = v->dist + 1;
 229                     }
 230                 }
 231                 if( e0 > 0 )
 232                     break;
 233             }
 234             // exclude the vertex from the active list
 235             first = first->next;
 236             v->next = 0;
 237         }
 238
 239         if( e0 <= 0 )
 240             break;
 241
 242         // find the minimum edge weight along the path
 243         minWeight = edgePtr[e0].weight;
 244         CV_Assert( minWeight > 0 );
 245         // k = 1: source tree, k = 0: destination tree
 246         for( int k = 1; k >= 0; k-- )
 247         {
 248             for( v = vtxPtr+edgePtr[e0^k].dst;; v = vtxPtr+edgePtr[ei].dst )
 249             {
 250                 if( (ei = v->parent) < 0 )
 251                     break;
 252                 weight = edgePtr[ei^k].weight;
 253                 minWeight = MIN(minWeight, weight);
 254                 CV_Assert( minWeight > 0 );
 255             }
 256             weight = fabs(v->weight);
 257             minWeight = MIN(minWeight, weight);
 258             CV_Assert( minWeight > 0 );
 259         }
 260
 261         // modify weights of the edges along the path and collect orphans
 262         edgePtr[e0].weight -= minWeight;
 263         edgePtr[e0^1].weight += minWeight;
 264         flow += minWeight;
 265
 266         // k = 1: source tree, k = 0: destination tree
 267         for( int k = 1; k >= 0; k-- )
 268         {
 269             for( v = vtxPtr+edgePtr[e0^k].dst;; v = vtxPtr+edgePtr[ei].dst )
 270             {
 271                 if( (ei = v->parent) < 0 )
 272                     break;
 273                 edgePtr[ei^(k^1)].weight += minWeight;
 274                 if( (edgePtr[ei^k].weight -= minWeight) == 0 )
 275                 {
 276                     orphans.push_back(v);
 277                     v->parent = ORPHAN;
 278                 }
 279             }
 280
 281             v->weight = v->weight + minWeight*(1-k*2);
 282             if( v->weight == 0 )
 283             {
 284                orphans.push_back(v);
 285                v->parent = ORPHAN;
 286             }
 287         }
 288
 289         // restore the search trees by finding new parents for the orphans
 290         curr_ts++;
 291         while( !orphans.empty() )
 292         {
 293             Vtx* v2 = orphans.back();
 294             orphans.pop_back();
 295
 296             int d, minDist = INT_MAX;
 297             e0 = 0;
 298             vt = v2->t;
 299
 300             for( ei = v2->first; ei != 0; ei = edgePtr[ei].next )
 301             {
 302                 if( edgePtr[ei^(vt^1)].weight == 0 )
 303                     continue;
 304                 u = vtxPtr+edgePtr[ei].dst;
 305                 if( u->t != vt || u->parent == 0 )
 306                     continue;
 307                 // compute the distance to the tree root
 308                 for( d = 0;; )
 309                 {
 310                     if( u->ts == curr_ts )
 311                     {
 312                         d += u->dist;
 313                         break;
 314                     }
 315                     ej = u->parent;
 316                     d++;
 317                     if( ej < 0 )
 318                     {
 319                         if( ej == ORPHAN )
 320                             d = INT_MAX-1;
 321                         else
 322                         {
 323                             u->ts = curr_ts;
 324                             u->dist = 1;
 325                         }
 326                         break;
 327                     }
 328                     u = vtxPtr+edgePtr[ej].dst;
 329                 }
 330
 331                 // update the distance
 332                 if( ++d < INT_MAX )
 333                 {
 334                     if( d < minDist )
 335                     {
 336                         minDist = d;
 337                         e0 = ei;
 338                     }
 339                     for( u = vtxPtr+edgePtr[ei].dst; u->ts != curr_ts; u = vtxPtr+edgePtr[u->parent].dst )
 340                     {
 341                         u->ts = curr_ts;
 342                         u->dist = --d;
 343                     }
 344                 }
 345             }
 346
 347             if( (v2->parent = e0) > 0 )
 348             {
 349                 v2->ts = curr_ts;
 350                 v2->dist = minDist;
 351                 continue;
 352             }
 353
 354             /* no parent is found */
 355             v2->ts = 0;
 356             for( ei = v2->first; ei != 0; ei = edgePtr[ei].next )
 357             {
 358                 u = vtxPtr+edgePtr[ei].dst;
 359                 ej = u->parent;
 360                 if( u->t != vt || !ej )
 361                     continue;
 362                 if( edgePtr[ei^(vt^1)].weight && !u->next )
 363                 {
 364                     u->next = nilNode;
 365                     last = last->next = u;
 366                 }
 367                 if( ej > 0 && vtxPtr+edgePtr[ej].dst == v2 )
 368                 {
 369                     orphans.push_back(u);
 370                     u->parent = ORPHAN;
 371                 }
 372             }
 373         }
 374     }
 375     return flow;
 376 }
 377
 378 template <class TWeight>
 379 bool GCGraph<TWeight>::inSourceSegment( int i )
 380 {
 381     CV_Assert( i>=0 && i<(int)vtcs.size() );
 382     return vtcs[i].t == 0;
 383 }
 384
 385 #endif