1/* 2 Maxflow package implements the max-flow(min-cuts, graph-cuts) algorithm that is used to solve the labeling problem in computer vision or graphics area. 3 4 The algorithm is described in 5 6 An Experimental Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Computer Vision. 7 Yuri Boykov and Vladimir Kolmogorov. 8 In IEEE Transactions on Pattern Analysis and Machine Intelligence, September 2004. 9 10 Reference the document of Graph struct for usage information. 11*/ 12package maxflow 13 14// import "fmt" 15 16type CapType int 17 18type arc struct { 19 head *Node /* node the arc points to */ 20 next *arc /* next arc with the same originating node */ 21 sister *arc /* reverse arc */ 22 rCap CapType /* residual capacity */ 23} 24 25type Node struct { 26 first *arc /* first outcoming arc */ 27 parent *arc /* node's parent */ 28 next *Node /* pointer to the next active node (or to itself if it is the last node in the list) */ 29 30 dist int /* distance to the terminal */ 31 counter int /* timestamp showing when dist was computed */ 32 33 isSink bool /* flag showing whether the node is in the source or in the sink tree */ 34 trCap CapType /* if tr_cap > 0 then tr_cap is residual capacity of the arc SOURCE->node 35 otherwise -tr_cap is residual capacity of the arc node->SINK */ 36} 37 38var terminalArc *arc = &arc{} 39var orphanArc *arc = &arc{} 40 41const INFINITE_D int = 1000000000 42 43type nodePtr struct { 44 ptr *Node 45 next *nodePtr 46} 47 48/* 49Graph is a data structure representing a graph for maxflow algorithm. 50 51Usage: 52 g := NewGraph() 53 54 nodes := make([]*Node, 2) 55 56 for i := range(nodes) { 57 nodes[i] = g.AddNode() 58 } // for i 59 60 g.SetTweights(nodes[0], 1, 5) 61 g.SetTweights(nodes[1], 2, 6) 62 g.AddEdge(nodes[0], nodes[1], 3, 4) 63 64 g.Run(); 65 66 flow := g.Flow() 67 68 if g.IsSource(nodes[0]) { 69 fmt.Println("nodes 0 is SOURCE") 70 } else { 71 fmt.Println("nodes 0 is SINK") 72 } // else 73*/ 74type Graph struct { 75 nodes []*Node 76 77 flow CapType /* total flow */ 78 queueFirst, queueLast *Node 79 orphanFirst, orphanLast *nodePtr 80 counter int 81 finish bool 82} 83 84// NewGraph creates an initialzed Graph instance. 85func NewGraph() *Graph { 86 return &Graph{} 87} 88 89// AddNode creates a node in the graph and returns a pointer to the node. 90// 91// Fields of the Node struct is not(and need not be) accessible. 92func (g *Graph) AddNode() *Node { 93 nd := &Node{} 94 g.nodes = append(g.nodes, nd) 95 return nd 96} 97 98// SetTweights sets the capacities of a node to the souce and sink node 99// 100// Do not call this method twice for a node 101func (g *Graph) SetTweights(i *Node, capSource, capSink CapType) { 102 if capSource < capSink { 103 g.flow += capSource 104 } else { 105 g.flow += capSink 106 } // else 107 108 i.trCap = capSource - capSink 109} 110 111// AddEdge adds edges between two nodes. 112// 113// cap and revCap are two directions 114func (g *Graph) AddEdge(from, to *Node, cap, revCap CapType) { 115 a, aRev := &arc{}, &arc{} 116 117 a.sister, aRev.sister = aRev, a 118 119 a.next = from.first 120 from.first = a 121 122 aRev.next = to.first 123 to.first = aRev 124 125 a.head, aRev.head = to, from 126 a.rCap, aRev.rCap = cap, revCap 127} 128 129// Flow returns the calculated maxflow. 130// 131// Call this method after calling to Run 132func (g *Graph) Flow() CapType { 133 return g.flow 134} 135 136// IsSource checks whether a node is a source in the minimum cuts. 137// 138// Call this method after calling to Run 139func (g *Graph) IsSource(i *Node) bool { 140 return i.parent != nil && !i.isSink 141} 142 143// Run executes the maxflow algorithm to find the maxflow of the current graph. 144func (g *Graph) Run() { 145 if g.finish { 146 return 147 } // if 148 149 var j, currentNode *Node 150 var a *arc 151 var np, npNext *nodePtr 152 153 g.maxflowInit() 154 155 for { 156 i := currentNode 157 if i != nil { 158 i.next = nil 159 if i.parent == nil { 160 i = nil 161 } // if 162 } // if 163 164 if i == nil { 165 i = g.nextActive() 166 if i == nil { 167 break 168 } // if 169 } // if 170 171 /* growth */ 172 if !i.isSink { 173 /* grow source tree */ 174 for a = i.first; a != nil; a = a.next { 175 if a.rCap != 0 { 176 j = a.head 177 if j.parent == nil { 178 j.isSink = false 179 j.parent = a.sister 180 j.counter = i.counter 181 j.dist = i.dist + 1 182 g.setActive(j) 183 } else if j.isSink { 184 break 185 } else if j.counter <= i.counter && j.dist > i.dist { 186 /* heuristic - trying to make the distance from j to the source shorter */ 187 j.parent = a.sister 188 j.counter = i.counter 189 j.dist = i.dist + 1 190 } // else if 191 } // if 192 } // for a 193 } else { 194 /* grow sink tree */ 195 for a = i.first; a != nil; a = a.next { 196 if a.sister.rCap != 0 { 197 j = a.head 198 if j.parent == nil { 199 j.isSink = true 200 j.parent = a.sister 201 j.counter = i.counter 202 j.dist = i.dist + 1 203 g.setActive(j) 204 } else if !j.isSink { 205 a = a.sister 206 break 207 } else if j.counter <= i.counter && j.dist > i.dist { 208 /* heuristic - trying to make the distance from j to the sink shorter */ 209 j.parent = a.sister 210 j.counter = i.counter 211 j.dist = i.dist + 1 212 } // else if 213 } // if 214 } // for a 215 } // else 216 217 g.counter++ 218 219 if a != nil { 220 /* set active flag */ 221 i.next = i 222 currentNode = i 223 224 /* augmentation */ 225 g.augment(a) 226 /* augmentation end */ 227 228 /* adoption */ 229 for np = g.orphanFirst; np != nil; np = g.orphanFirst { 230 npNext = np.next 231 np.next = nil 232 233 for np = g.orphanFirst; np != nil; np = g.orphanFirst { 234 //nodeptr_block -> Delete(np); TODO reuse nodeptr 235 g.orphanFirst = np.next 236 i = np.ptr 237 if g.orphanFirst == nil { 238 g.orphanLast = nil 239 } // if 240 if i.isSink { 241 g.processSinkOrphan(i) 242 } else { 243 g.processSourceOrphan(i) 244 } // else 245 } // for np 246 247 g.orphanFirst = npNext 248 } // for np 249 /* adoption end */ 250 } else { 251 currentNode = nil 252 } // else 253 } // for true 254 255 g.finish = true 256} 257 258func (g *Graph) maxflowInit() { 259 g.queueFirst, g.queueLast = nil, nil 260 261 for _, i := range g.nodes { 262 i.next = nil 263 i.counter = 0 264 if i.trCap > 0 { 265 /* i is connected to the source */ 266 i.isSink = false 267 i.parent = terminalArc 268 g.setActive(i) 269 i.counter = 0 270 i.dist = 1 271 } else if i.trCap < 0 { 272 /* i is connected to the sink */ 273 i.isSink = true 274 i.parent = terminalArc 275 g.setActive(i) 276 i.counter = 0 277 i.dist = 1 278 } else { 279 i.parent = nil 280 } // else 281 } // for i 282 g.counter = 0 283} 284 285/* 286 nextActive returns the next active node. 287 If it is connected to the sink, it stays in the list, 288 otherwise it is removed from the list 289*/ 290func (g *Graph) nextActive() *Node { 291 for { 292 i := g.queueFirst 293 if i == nil { 294 return nil 295 } // if 296 297 /* remove it from the active list */ 298 if i.next == i { 299 g.queueFirst, g.queueLast = nil, nil 300 } else { 301 g.queueFirst = i.next 302 } // else 303 i.next = nil 304 305 /* a node in the list is active iff it has a parent */ 306 if i.parent != nil { 307 return i 308 } // if 309 } // for true 310 311 return nil 312} 313 314/* 315 Functions for processing active list. 316 i->next points to the next node in the list (or to i, if i is the last node in the list). 317 If i->next is NULL iff i is not in the list. 318*/ 319func (g *Graph) setActive(i *Node) { 320 if i.next == nil { 321 /* it's not in the list yet */ 322 if g.queueLast != nil { 323 g.queueLast.next = i 324 } else { 325 g.queueFirst = i 326 } // else 327 g.queueLast = i 328 i.next = i 329 } // if 330} 331 332func (g *Graph) augment(middleArc *arc) { 333 var i *Node 334 var a *arc 335 var bottleNeck CapType 336 var np *nodePtr 337 338 /* 1. Finding bottleneck capacity */ 339 /* 1a - the source tree */ 340 bottleNeck = middleArc.rCap 341 for i = middleArc.sister.head; true; i = a.head { 342 a = i.parent 343 if a == terminalArc { 344 break 345 } // if 346 if bottleNeck > a.sister.rCap { 347 bottleNeck = a.sister.rCap 348 } // if 349 } // for i 350 if bottleNeck > i.trCap { 351 bottleNeck = i.trCap 352 } // if 353 /* 1b - the sink tree */ 354 for i = middleArc.head; true; i = a.head { 355 a = i.parent 356 if a == terminalArc { 357 break 358 } // if 359 if bottleNeck > a.rCap { 360 bottleNeck = a.rCap 361 } // if 362 } // for i 363 if bottleNeck > -i.trCap { 364 bottleNeck = -i.trCap 365 } // if 366 367 /* 2. Augmenting */ 368 /* 2a - the source tree */ 369 middleArc.sister.rCap += bottleNeck 370 middleArc.rCap -= bottleNeck 371 for i = middleArc.sister.head; true; i = a.head { 372 a = i.parent 373 if a == terminalArc { 374 break 375 } // if 376 a.rCap += bottleNeck 377 a.sister.rCap -= bottleNeck 378 if a.sister.rCap == 0 { 379 /* add i to the adoption list */ 380 i.parent = orphanArc 381 np = &nodePtr{} 382 np.ptr = i 383 np.next = g.orphanFirst 384 g.orphanFirst = np 385 } // if 386 } // for i 387 i.trCap -= bottleNeck 388 if i.trCap == 0 { 389 /* add i to the adoption list */ 390 i.parent = orphanArc 391 np = &nodePtr{} 392 np.ptr = i 393 np.next = g.orphanFirst 394 g.orphanFirst = np 395 } // if 396 /* 2b - the sink tree */ 397 for i = middleArc.head; true; i = a.head { 398 a = i.parent 399 if a == terminalArc { 400 break 401 } // if 402 a.sister.rCap += bottleNeck 403 a.rCap -= bottleNeck 404 if a.rCap == 0 { 405 /* add i to the adoption list */ 406 i.parent = orphanArc 407 np = &nodePtr{} 408 np.ptr = i 409 np.next = g.orphanFirst 410 g.orphanFirst = np 411 } // if 412 } // for i 413 i.trCap += bottleNeck 414 if i.trCap == 0 { 415 /* add i to the adoption list */ 416 i.parent = orphanArc 417 np = &nodePtr{} 418 np.ptr = i 419 np.next = g.orphanFirst 420 g.orphanFirst = np 421 } // if 422 423 g.flow += bottleNeck 424} 425 426func (g *Graph) processSinkOrphan(i *Node) { 427 var a0Min *arc 428 var dMin int = INFINITE_D 429 430 /* trying to find a new parent */ 431 for a0 := i.first; a0 != nil; a0 = a0.next { 432 if a0.rCap != 0 { 433 j := a0.head 434 if a := j.parent; j.isSink && a != nil { 435 /* checking the origin of j */ 436 //d = 0; 437 var d int = 0 438 for true { 439 if j.counter == g.counter { 440 d += j.dist 441 break 442 } // if 443 a = j.parent 444 d++ 445 if a == terminalArc { 446 j.counter = g.counter 447 j.dist = 1 448 break 449 } // if 450 if a == orphanArc { 451 d = INFINITE_D 452 break 453 } // if 454 j = a.head 455 } // for true 456 /* j originates from the sink - done */ 457 if d < INFINITE_D { 458 if d < dMin { 459 a0Min = a0 460 dMin = d 461 } // if 462 /* set marks along the path */ 463 for j := a0.head; j.counter != g.counter; j = j.parent.head { 464 j.counter = g.counter 465 j.dist = d 466 d-- 467 } // for j 468 } // if 469 } // if 470 } // if 471 } // for a0 472 473 if i.parent = a0Min; i.parent != nil { 474 i.counter = g.counter 475 i.dist = dMin + 1 476 } else { 477 /* no parent is found */ 478 i.counter = 0 479 480 /* process neighbors */ 481 for a0 := i.first; a0 != nil; a0 = a0.next { 482 j := a0.head 483 if a := j.parent; j.isSink && a != nil { 484 if a0.rCap != 0 { 485 g.setActive(j) 486 } // if 487 if a != terminalArc && a != orphanArc && a.head == i { 488 /* add j to the adoption list */ 489 i.parent = orphanArc 490 np := &nodePtr{} 491 np.ptr = j 492 if g.orphanLast != nil { 493 g.orphanLast.next = np 494 } else { 495 g.orphanFirst = np 496 } // else 497 g.orphanLast = np 498 np.next = nil 499 } // i f 500 } // if 501 } // for a0 502 } // else 503} 504 505func (g *Graph) processSourceOrphan(i *Node) { 506 var a0Min *arc 507 var dMin int = INFINITE_D 508 509 /* trying to find a new parent */ 510 for a0 := i.first; a0 != nil; a0 = a0.next { 511 if a0.sister.rCap != 0 { 512 j := a0.head 513 if a := j.parent; j.isSink && a != nil { 514 /* checking the origin of j */ 515 var d int = 0 516 for true { 517 if j.counter == g.counter { 518 d += j.dist 519 break 520 } // if 521 a = j.parent 522 d++ 523 if a == terminalArc { 524 j.counter = g.counter 525 j.dist = 1 526 break 527 } // if 528 if a == orphanArc { 529 d = INFINITE_D 530 break 531 } // if 532 j = a.head 533 } // for true 534 /* j originates from the source - done */ 535 if d < INFINITE_D { 536 if d < dMin { 537 a0Min = a0 538 dMin = d 539 } // if 540 /* set marks along the path */ 541 for j := a0.head; j.counter != g.counter; j = j.parent.head { 542 j.counter = g.counter 543 j.dist = d 544 d-- 545 } // for j 546 } // if 547 } // if 548 } // if 549 } // for a0 550 551 if i.parent = a0Min; i.parent != nil { 552 i.counter = g.counter 553 i.dist = dMin + 1 554 } else { 555 /* no parent is found */ 556 i.counter = 0 557 558 /* process neighbors */ 559 for a0 := i.first; a0 != nil; a0 = a0.next { 560 j := a0.head 561 if a := j.parent; !j.isSink && a != nil { 562 if a0.sister.rCap != 0 { 563 g.setActive(j) 564 } // if 565 if a != terminalArc && a != orphanArc && a.head == i { 566 /* add j to the adoption list */ 567 j.parent = orphanArc 568 np := &nodePtr{} 569 np.ptr = j 570 if g.orphanLast != nil { 571 g.orphanLast.next = np 572 } else { 573 g.orphanFirst = np 574 } // else 575 g.orphanLast = np 576 np.next = nil 577 } // if 578 } // if 579 } // for a0 580 } // else 581} 582