Graph cut optimization

Short description: Combinatorial optimization method for a family of functions of discrete variables

Graph cut optimization is a combinatorial optimization method applicable to a family of functions of discrete variables, named after the concept of cut in the theory of flow networks. Thanks to the max-flow min-cut theorem, determining the minimum cut over a graph representing a flow network is equivalent to computing the maximum flow over the network. Given a pseudo-Boolean function [math]\displaystyle{ f }[/math], if it is possible to construct a flow network with positive weights such that

• each cut [math]\displaystyle{ C }[/math] of the network can be mapped to an assignment of variables [math]\displaystyle{ \mathbf{x} }[/math] to [math]\displaystyle{ f }[/math] (and vice versa), and
• the cost of [math]\displaystyle{ C }[/math] equals [math]\displaystyle{ f(\mathbf{x}) }[/math] (up to an additive constant)

then it is possible to find the global optimum of [math]\displaystyle{ f }[/math] in polynomial time by computing a minimum cut of the graph. The mapping between cuts and variable assignments is done by representing each variable with one node in the graph and, given a cut, each variable will have a value of 0 if the corresponding node belongs to the component connected to the source, or 1 if it belong to the component connected to the sink.

Not all pseudo-Boolean functions can be represented by a flow network, and in the general case the global optimization problem is NP-hard. There exist sufficient conditions to characterise families of functions that can be optimised through graph cuts, such as submodular quadratic functions. Graph cut optimization can be extended to functions of discrete variables with a finite number of values, that can be approached with iterative algorithms with strong optimality properties, computing one graph cut at each iteration.

Graph cut optimization is an important tool for inference over graphical models such as Markov random fields or conditional random fields, and it has applications in computer vision problems such as image segmentation, denoising, registration and stereo matching.

Representability

A pseudo-Boolean function [math]\displaystyle{ f: \{0, 1\}^n \to \mathbb{R} }[/math] is said to be representable if there exists a graph [math]\displaystyle{ G = (V, E) }[/math] with non-negative weights and with source and sink nodes [math]\displaystyle{ s }[/math] and [math]\displaystyle{ t }[/math] respectively, and there exists a set of nodes [math]\displaystyle{ V_0 = \{v_1, \dots, v_n\} \subset V - \{s, t\} }[/math] such that, for each tuple of values [math]\displaystyle{ (x_1, \dots, x_n) \in \{0, 1\}^n }[/math] assigned to the variables, [math]\displaystyle{ f(x_1, \dots, x_n) }[/math] equals (up to a constant) the value of the flow determined by a minimum cut [math]\displaystyle{ C = (S, T) }[/math] of the graph [math]\displaystyle{ G }[/math] such that [math]\displaystyle{ v_i \in S }[/math] if [math]\displaystyle{ x_i = 0 }[/math] and [math]\displaystyle{ v_i \in T }[/math] if [math]\displaystyle{ x_i = 1 }[/math].

It is possible to classify pseudo-Boolean functions according to their order, determined by the maximum number of variables contributing to each single term. All first order functions, where each term depends upon at most one variable, are always representable. Quadratic functions

[math]\displaystyle{ f(\mathbf{x}) = w_0 + \sum_i w_i(x_i) + \sum_{i \lt j} w_{ij}(x_i, x_j) . }[/math]

are representable if and only if they are submodular, i.e. for each quadratic term [math]\displaystyle{ w_{ij} }[/math] the following condition is satisfied

[math]\displaystyle{ w_{ij}(0, 0) + w_{ij}(1, 1) \le w_{ij}(0, 1) + w_{ij}(1, 0) . }[/math]

Cubic functions

[math]\displaystyle{ f(\mathbf{x}) = w_0 + \sum_i w_i(x_i) + \sum_{i \lt j} w_{ij}(x_i, x_j) + \sum_{i \lt j \lt k} w_{ijk}(x_i, x_j, x_k) }[/math]

are representable if and only if they are regular, i.e. all possible binary projections to two variables, obtained by fixing the value of the remaining variable, are submodular. For higher-order functions, regularity is a necessary condition for representability.

Graph construction

Graph construction for a representable function is simplified by the fact that the sum of two representable functions [math]\displaystyle{ f' }[/math] and [math]\displaystyle{ f'' }[/math] is representable, and its graph [math]\displaystyle{ G = (V' \cup V'', E' \cup E'') }[/math] is the union of the graphs [math]\displaystyle{ G' = (V', E') }[/math] and [math]\displaystyle{ G'' = (V'', E'') }[/math] representing the two functions. Such theorem allows to build separate graphs representing each term and combine them to obtain a graph representing the entire function.

The graph representing a quadratic function of [math]\displaystyle{ n }[/math] variables contains [math]\displaystyle{ n + 2 }[/math] vertices, two of them representing the source and sink and the others representing the variables. When representing higher-order functions, the graph contains auxiliary nodes that allow to model higher-order interactions.

Unary terms

A unary term [math]\displaystyle{ w_i }[/math] depends only on one variable [math]\displaystyle{ x_i }[/math] and can be represented by a graph with one non-terminal node [math]\displaystyle{ v_i }[/math] and one edge [math]\displaystyle{ s \rightarrow v_i }[/math] with weight [math]\displaystyle{ w_i(1) - w_i(0) }[/math] if [math]\displaystyle{ w_i(1) \ge w_i(0) }[/math], or [math]\displaystyle{ v_i \rightarrow t }[/math] with weight [math]\displaystyle{ w_i(0) - w_i(1) }[/math] if [math]\displaystyle{ w_i(1) \lt w_i(0) }[/math].

Binary terms

Example of a graph representing a quadratic term [math]\displaystyle{ w_{ij}(x_i, x_j) }[/math] in case [math]\displaystyle{ w_{ij}(1, 0) - w_{ij}(0, 0) \gt 0 }[/math] and [math]\displaystyle{ w_{ij}(1, 1) - w_{ij}(1, 0) \lt 0 }[/math].

A quadratic (or binary) term [math]\displaystyle{ w_{ij} }[/math] can be represented by a graph containing two non-terminal nodes [math]\displaystyle{ v_i }[/math] and [math]\displaystyle{ v_j }[/math]. The term can be rewritten as

[math]\displaystyle{ w_{ij}(x_i, x_j) = w_{ij}(0, 0) + k_i x_i + k_j x_j + k_{ij} \left( (1 - x_i) x_j + x_i (1 - x_j) \right) }[/math]

with

[math]\displaystyle{ \begin{align} k_i &= \frac{1}{2} (w_{ij}(1, 0) - w_{ij}(0, 0)) \\ k_j &= \frac{1}{2} (w_{ij}(1, 1) - w_{ij}(1, 0)) \\ k_{ij} &= \frac{1}{2} (w_{ij}(0, 1) + w_{ij}(1, 0) - w_{ij}(0, 0) - w_{ij}(1, 1)) . \end{align} }[/math]

In this expression, the first term is constant and it is not represented by any edge, the two following terms depend on one variable and are represented by one edge, as shown in the previous section for unary terms, while the third term is represented by an edge [math]\displaystyle{ v_i \rightarrow v_j }[/math] with weight [math]\displaystyle{ w_{ij}(0, 1) + w_{ij}(1, 0) - w_{ij}(0, 0) - w_{ij}(1, 1) }[/math] (submodularity guarantees that the weight is non-negative).

Ternary terms

A cubic (or ternary) term [math]\displaystyle{ w_{ijk} }[/math] can be represented by a graph with four non-terminal nodes, three of them ([math]\displaystyle{ v_i }[/math], [math]\displaystyle{ v_j }[/math] and [math]\displaystyle{ v_k }[/math]) associated to the three variables plus one fourth auxiliary node [math]\displaystyle{ v_{ijk} }[/math]. A generic ternary term can be rewritten as the sum of a constant, three unary terms, three binary terms, and a ternary term in simplified form. There may be two different cases, according to the sign of [math]\displaystyle{ p = w_{ijk}(0, 0, 0) + w_{ijk}(0, 1, 1) + w_{ijk}(1, 0, 1) + w_{ijk}(1, 1, 0) }[/math]. If [math]\displaystyle{ p \gt 0 }[/math] then

[math]\displaystyle{ w_{ijk}(x_i, x_j, x_k) = w_{ijk}(0, 0, 0) + p_1 (x_i - 1) + p_2 (x_j - 1) + p_3 (x_k - 1) + p_{23}(x_j - 1) x_k + p_{31} x_i (x_k - 1) + p_{12} (x_i - 1) x_j - p x_i x_j x_k }[/math]

Example of a graph representing the ternary term [math]\displaystyle{ p x_i x_j x_k }[/math] when [math]\displaystyle{ p \gt 0 }[/math] (left) and when [math]\displaystyle{ p \lt 0 }[/math] (right).

with

[math]\displaystyle{ \begin{align} p_1 &= w_{ijk}(1, 0, 1) - w_{ijk}(0, 0, 1) \\ p_2 &= w_{ijk}(1, 1, 0) - w_{ijk}(1, 0, 1) \\ p_3 &= w_{ijk}(0, 1, 1) - w_{ijk}(0, 1, 0) \\ p_{23} &= w_{ijk}(0, 0, 1) + w_{ijk}(0, 1, 0) - w_{ijk}(0, 0, 0) - w_{ijk}(0, 1, 1) \\ p_{31} &= w_{ijk}(0, 0, 1) + w_{ijk}(1, 0, 0) - w_{ijk}(0, 0, 0) - w_{ijk}(1, 0, 1) \\ p_{12} &= w_{ijk}(0, 1, 0) + w_{ijk}(1, 0, 0) - w_{ijk}(0, 0, 0) - w_{ijk}(1, 1, 0) . \end{align} }[/math]

If [math]\displaystyle{ p \lt 0 }[/math] the construction is similarly, but the variables will have opposite value. If the function is regular, then all its projections of two variables will be submodular, implying that [math]\displaystyle{ p_{23} }[/math], [math]\displaystyle{ p_{31} }[/math] and [math]\displaystyle{ p_{12} }[/math] are positive and then all terms in the new representation are submodular.

In this decomposition, the constant, unary and binary terms can be represented as shown in the previous sections. If [math]\displaystyle{ p \gt 0 }[/math] the ternary term can be represented with a graph with four edges [math]\displaystyle{ v_i \rightarrow v_{ijk} }[/math], [math]\displaystyle{ v_j \rightarrow v_{ijk} }[/math], [math]\displaystyle{ v_k \rightarrow v_{ijk} }[/math], [math]\displaystyle{ v_{ijk} \rightarrow t }[/math], all with weight [math]\displaystyle{ p }[/math], while if [math]\displaystyle{ p \lt 0 }[/math] the term can be represented by four edges [math]\displaystyle{ v_{ijk} \rightarrow v_i }[/math], [math]\displaystyle{ v_{ijk} \rightarrow v_j }[/math], [math]\displaystyle{ v_{ijk} \rightarrow v_k }[/math], [math]\displaystyle{ s \rightarrow v_{ijk} }[/math] with weight [math]\displaystyle{ -p }[/math].

Minimum cut

After building a graph representing a pseudo-Boolean function, it is possible to compute a minimum cut using one among the various algorithms developed for flow networks, such as Ford–Fulkerson, Edmonds–Karp, and Boykov–Kolmogorov algorithm. The result is a partition of the graph in two connected components [math]\displaystyle{ S }[/math] and [math]\displaystyle{ T }[/math] such that [math]\displaystyle{ s \in S }[/math] and [math]\displaystyle{ t \in T }[/math], and the function attains its global minimum when [math]\displaystyle{ x_i = 0 }[/math] for each [math]\displaystyle{ i }[/math] such that the corresponding node [math]\displaystyle{ v_i \in S }[/math], and [math]\displaystyle{ x_i = 1 }[/math] for each [math]\displaystyle{ i }[/math] such that the corresponding node [math]\displaystyle{ v_i \in T }[/math].

Max-flow algorithms such as Boykov–Kolmogorov's are very efficient in practice for sequential computation, but they are difficult to parallelise, making them not suitable for distributed computing applications and preventing them from exploiting the potential of modern CPUs. Parallel max-flow algorithms were developed, such as push-relabel and jump-flood, that can also take advantage of hardware acceleration in GPGPU implementations.

Functions of discrete variables with more than two values

The previous construction allows global optimization of pseudo-Boolean functions only, but it can be extended to quadratic functions of discrete variables with a finite number of values, in the form

[math]\displaystyle{ f(\mathbf{x}) = \sum_{i \in V} D(x_i) + \sum_{(i, j) \in E} S(x_i, x_j) }[/math]

where [math]\displaystyle{ E \subseteq V \times V }[/math] and [math]\displaystyle{ x_i \in \Lambda = \{1, \dots, k\} }[/math]. The function [math]\displaystyle{ D(x_i) }[/math] represents the unary contribution of each variable (often referred as data term), while the function [math]\displaystyle{ S(x_i, x_j) }[/math] represents binary interactions between variables (smoothness term). In the general case, optimization of such functions is a NP-hard problem, and stochastic optimization methods such as simulated annealing are sensitive to local minima and in practice they can generate arbitrarily sub-optimal results. With graph cuts it is possible to construct move-making algorithms that allow to reach in polynomial time a local minima with strong optimality properties for a wide family of quadratic functions of practical interest (when the binary interaction [math]\displaystyle{ S(x_i, x_j) }[/math] is a metric or a semimetric), such that the value of the function in the solution lies within a constant and known factor from the global optimum.

Given a function [math]\displaystyle{ f: \Lambda^n \to \mathbb{R} }[/math] with [math]\displaystyle{ \Lambda = \{1, \dots, k\} }[/math], and a certain assignment of values [math]\displaystyle{ \mathbf{x} = (x_1, \dots, x_n) \in \Lambda^n }[/math] to the variables, it is possible to associate each assignment [math]\displaystyle{ \mathbf{x} }[/math] to a partition [math]\displaystyle{ P = \{P_l | l \in \Lambda \} }[/math] of the set of variables, such that, [math]\displaystyle{ P_l = \{ x_i | x_i = l \in \Lambda \} }[/math]. Give two distinct assignments [math]\displaystyle{ P }[/math] and [math]\displaystyle{ P' }[/math] and a value [math]\displaystyle{ \alpha \in \Lambda }[/math], a move that transforms [math]\displaystyle{ P }[/math] into [math]\displaystyle{ P' }[/math] is said to be an [math]\displaystyle{ \alpha }[/math]-expansion if [math]\displaystyle{ P_\alpha \subset P'_\alpha }[/math] and [math]\displaystyle{ P'_l \subset P_l \; \forall l \in \Lambda - \{ \alpha \} }[/math]. Given a couple of values [math]\displaystyle{ \alpha }[/math] and [math]\displaystyle{ \beta }[/math], a move is said to be an [math]\displaystyle{ \alpha\beta }[/math]-swap if [math]\displaystyle{ P_l = P'_l \; \forall l \in \Lambda - \{ \alpha, \beta \} }[/math]. Intuitively, an [math]\displaystyle{ \alpha }[/math]-expansion move from [math]\displaystyle{ \mathbf{x} }[/math] assigns the value of [math]\displaystyle{ \alpha }[/math] to some variables that have a different value in [math]\displaystyle{ \mathbf{x} }[/math], while an [math]\displaystyle{ \alpha\beta }[/math]-swap move assigns [math]\displaystyle{ \alpha }[/math] to some variables that have value [math]\displaystyle{ \beta }[/math] in [math]\displaystyle{ \mathbf{x} }[/math] and vice versa.

For each iteration, the [math]\displaystyle{ \alpha }[/math]-expansion algorithm computes, for each possible value [math]\displaystyle{ \alpha }[/math], the minimum of the function among all assignments [math]\displaystyle{ \Alpha(\mathbf{x}) }[/math] that can be reached with a single [math]\displaystyle{ \alpha }[/math]-expansion move from the current temporary solution [math]\displaystyle{ \mathbf{x} }[/math], and takes it as the new temporary solution.

[math]\displaystyle{ \mathbf{x} := \text{arbitrary value in } \Lambda^n }[/math]
[math]\displaystyle{ \text{exit} := 0 }[/math]
while [math]\displaystyle{ \text{exit} \ne 1 }[/math]:
    [math]\displaystyle{ \text{exit} = 1 }[/math]
    foreach [math]\displaystyle{ \alpha \in \Lambda }[/math]:
        [math]\displaystyle{ \mathbf{\hat{x}} := \arg \min_{\mathbf{y} \in \Alpha(\mathbf{x})} f(\mathbf{y}) }[/math]
        if [math]\displaystyle{ f(\mathbf{\hat{x}}) \lt  f(\mathbf{x}) }[/math]:
            [math]\displaystyle{ \mathbf{x} = \mathbf{\hat{x}} }[/math]
            [math]\displaystyle{ \text{exit} := 0 }[/math]

The [math]\displaystyle{ \alpha\beta }[/math]-swap algorithm is similar, but it searches for the minimum among all assignments [math]\displaystyle{ \Alpha\Beta(\mathbf{x}) }[/math] reachable with a single [math]\displaystyle{ \alpha\beta }[/math]-swap move from [math]\displaystyle{ \mathbf{x} }[/math].

[math]\displaystyle{ \mathbf{x} := \text{arbitrary value in } \Lambda^n }[/math]
[math]\displaystyle{ \text{exit} := 0 }[/math]
while [math]\displaystyle{ \text{exit} \ne 1 }[/math]:
    [math]\displaystyle{ \text{exit} = 1 }[/math]
    foreach [math]\displaystyle{ (\alpha, \beta) \in \Lambda^2 }[/math]:
        [math]\displaystyle{ \mathbf{\hat{x}} := \arg \min_{\mathbf{y} \in \Alpha\Beta(\mathbf{x})} f(\mathbf{y}) }[/math]
        if [math]\displaystyle{ f(\mathbf{\hat{x}}) \lt  f(\mathbf{x}) }[/math]:
            [math]\displaystyle{ \mathbf{x} = \mathbf{\hat{x}} }[/math]
            [math]\displaystyle{ \text{exit} := 0 }[/math]

In both cases, the optimization problem in the innermost loop can be solved exactly and efficiently with a graph cut. Both algorithms terminate certainly in a finite number of iterations of the outer loop, and in practice such number is small, with most of the improvement happening at the first iteration. The algorithms can generate different solutions depending on the initial guess, but in practice they are robust with respect to initialisation, and starting with a point where all variables are assigned to the same random value is usually sufficient to produce good quality results.

The solution generated by such algorithms is not necessarily a global optimum, but it has strong guarantees of optimality. If [math]\displaystyle{ S(x_i, x_j) }[/math] is a metric and [math]\displaystyle{ \mathbf{x} }[/math] is a solution generated by the [math]\displaystyle{ \alpha }[/math]-expansion algorithm, or if [math]\displaystyle{ S(x_i, x_j) }[/math] is a semimetric and [math]\displaystyle{ \mathbf{x} }[/math] is a solution generated by the [math]\displaystyle{ \alpha\beta }[/math]-swap algorithm, then [math]\displaystyle{ f(\mathbf{x}) }[/math] lies within a known and constant factor from the global minimum [math]\displaystyle{ f(\mathbf{x}^*) }[/math]:

[math]\displaystyle{ f(\mathbf{x}) \le 2 \frac{ \max_{\alpha \ne \beta \in \Lambda} S(\alpha, \beta) }{ \min_{\alpha \ne \beta \in \Lambda} S(\alpha, \beta) } f(\mathbf{x}^*) . }[/math]

Non-submodular functions

Generally speaking, the problem of optimizing a non-submodular pseudo-Boolean function is NP-hard and cannot be solved in polynomial time with a simple graph cut. The simplest approach is to approximate the function with a similar but submodular one, for instance truncating all non-submodular terms or replacing them with similar submodular expressions. Such approach is generally sub-optimal, and it produces acceptable results only if the number of non-submodular terms is relatively small.

In case of quadratic non-submodular functions, it is possible to compute in polynomial time a partial solution using algorithms such as QPBO. Higher-order functions can be reduced in polynomial time to a quadratic form that can be optimised with QPBO.

Higher-order functions

Quadratic functions are extensively studied and were characterised in detail, but more general results were derived also for higher-order functions. While quadratic functions can indeed model many problems of practical interest, they are limited by the fact they can represent only binary interactions between variables. The possibility to capture higher-order interactions allows to better capture the nature of the problem and it can provide higher quality results that could be difficult to achieve with quadratic models. For instance in computer vision applications, where each variable represents a pixel or voxel of the image, higher-order interactions can be used to model texture information, that would be difficult to capture using only quadratic functions.

Sufficient conditions analogous to submodularity were developed to characterise higher-order pseudo-Boolean functions that can be optimised in polynomial time, and there exists algorithms analogous to [math]\displaystyle{ \alpha }[/math]-expansion and [math]\displaystyle{ \alpha\beta }[/math]-swap for some families of higher-order functions. The problem is NP-hard in the general case, and approximate methods were developed for fast optimization of functions that do not satisfy such conditions.

References

• Boykov, Yuri; Veksler, Olga; Zabih, Ramin (2001). "Fast approximate energy minimization via graph cuts". IEEE Transactions on Pattern Analysis and Machine Intelligence 23 (11): 1222–1239. doi:10.1109/34.969114. 
• Freedman, Daniel; Drineas, Petros (2005). "Energy minimization via graph cuts: Settling what is possible". IEEE Computer Society Conference on Computer Vision and Pattern Recognition. 2. pp. 939–946. http://www.vision.cs.rpi.edu/publications/pdfs/freedman_cvpr05-cuts.pdf. 
• Goldberg, Andrew V; Tarjan, Robert E (1988). "A new approach to the maximum–flow problem". Journal of the ACM 35 (4): 921–940. doi:10.1145/48014.61051. http://akira.ruc.dk/~keld/teaching/algoritmedesign_f03/Artikler/08/Goldberg88.pdf. 
• Ishikawa, Hiroshi (2014). "Higher–Order Clique Reduction Without Auxiliary Variables". IEEE Conference on Computer Vision and Pattern Recognition. IEEE. pp. 1362–1369. https://www.cv-foundation.org/openaccess/content_cvpr_2014/papers/Ishikawa_Higher-Order_Clique_Reduction_2014_CVPR_paper.pdf. 
• Hong, Li; Chen, George (2004). "Segment-based stereo matching using graph cuts". Proceedings of the 2004 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. 1. pp. 74–81. http://vc.cs.nthu.edu.tw/home/paper/codfiles/qyliu/201211290607/Segment-based%20Stereo%20Matching%20Using%20Graph%20Cuts.pdf. 
• Kohli, Pushmeet; Kumar, M. Pawan; Torr, Philip H.S. (2009). "P³ & Beyond: Move Making Algorithms for Solving Higher Order Functions". IEEE Transactions on Pattern Analysis and Machine Intelligence 31 (9): 1645–1656. doi:10.1109/tpami.2008.217. PMID 19574624. http://www.robots.ox.ac.uk/~phst/Papers/2008/P3%20Journal/p3_pami.pdf. 
• Kolmogorov, Vladimir; Rother, Carsten (2007). "Minimizing Nonsubmodular Functions: A Review". IEEE Transactions on Pattern Analysis and Machine Intelligence 29 (7): 1274–1279. doi:10.1109/tpami.2007.1031. PMID 17496384. 
• Kolmogorov, Vladimir; Zabin, Ramin (2004). "What energy functions can be minimized via graph cuts?". IEEE Transactions on Pattern Analysis and Machine Intelligence 26 (2): 1645–1656. doi:10.1109/TPAMI.2004.1262177. PMID 15376891. https://ecommons.cornell.edu/bitstream/handle/1813/5842/2001-1857.pdf?sequence=1. 
• Lombaert, Herve; Cheriet, Farida (2012). "Simultaneous image de-noising and registration using graph cuts: Application to corrupted medical images". 11th International Conference on Information Science, Signal Processing and their Applications. pp. 264–268. http://cim.mcgill.ca/~lombaert/graphcuts-denoising-registration.pdf. 
• Peng, Yi; Chen, Li; Ou-Yang, Fang-Xin; Chen, Wei; Yong, Jun-Hai (2015). "JF-Cut: a parallel graph cut approach for large-scale image and video". IEEE Transactions on Image Processing 24 (2): 655–666. doi:10.1109/TIP.2014.2378060. PMID 25494510. Bibcode: 2015ITIP...24..655P. 
• Rother, Carsten; Kolmogorov, Vladimir; Blake, Andrew (2004). "Grabcut: Interactive foreground extraction using iterated graph cuts". ACM transactions on graphics. 23. pp. 309–314. http://pages.cs.wisc.edu/~dyer/cs534-fall11/papers/grabcut-rother.pdf. 
• So, Ronald WK; Tang, Tommy WH; Chung, Albert CS (2011). "Non-rigid image registration of brain magnetic resonance images using graph-cuts". Pattern Recognition 44 (10–11): 2450–2467. doi:10.1016/j.patcog.2011.04.008. 
• Stich, Timo (2009). "Graph Cuts with CUDA". GPU Technology Conference. https://www.nvidia.com/content/gtc/documents/1060_gtc09.pdf. 
• Tang, Tommy WH; Chung, Albert CS (2007). "Non-rigid image registration using graph-cuts". International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 916–924. doi:10.1007/978-3-540-75757-3_111. https://link.springer.com/content/pdf/10.1007/978-3-540-75757-3_111.pdf. 
• Vineet, Vibhav; Narayanan, PJ (2008). "CUDA cuts: Fast graph cuts on the GPU". IEEE Computer Society Conference on Computer Vision and Pattern Recognition Workshops. pp. 1–8. http://web2py.iiit.ac.in/publications/default/download/inproceedings.pdf.ac87b6f5-826a-41d1-91cf-b19f2574ac66.pdf. 

External links

• Implementation (C++) of several graph cut algorithms by Vladimir Kolmogorov.
• GCO, graph cut optimization library by Olga Veksler and Andrew Delong.