System and method for fast approximate energy minimization via graph cuts

[0002] This application claims priority of U.S. provisional applications Serial No. 60/151,413 entitled, “Fast Approximate Energy Minimization Via Graph Cuts” filed Aug. 30, 1999 by the present applicants.

[0004] This invention relates generally to computer vision and more particularly to minimizing energy functions in labeling pixels in early vision.

[0005] In computer vision, early vision is usually considered to involve the description of geometric structure in an image or sequence of images. The behavioral constraints on real-time visual systems typically require that the early vision stage be fast, reliable, general and automatic. Many early vision problems require estimating some spatially varying quantity such as intensity or disparity, from noisy measurements.

[0006] Spatially varying quantities tend to be piecewise smooth, i.e. they vary smoothly at most points, but change dramatically at object boundaries. Every pixel in a set P must be assigned a label in some set L. For motion or stereo, the labels are disparities. For image restoration, the labels represent intensities. The goal is to find a labeling f that assigns each pixel p∈P a label f_p∈L, where f is both piecewise smooth and consistent with observed data.

[0007] Computer vision problems can be formulated in terms of minimization of energy functions. Energy functions, however, are generally difficult to minimize. The major difficulty with energy minimization for computer vision lies in the enormous computational costs. Energy functions for computer vision typically have many local minima. Also, the space of possible labelings has the dimension |P| which is many thousands.

[0008] There have been many attempts to design fast algorithms for energy minimization. Due to the inefficiency of computing a global minimum, some solutions are directed instead at computing a local minimum. In general, a local minimum can be arbitrarily far from the optimum. It thus may not convey any of the global image properties that are encoded in the energy function. It is, however, difficult to determine the exact cause of an algorithmi's failures. When an algorithm gives unsatisfactory results, it may be due to either a poor choice of the energy function, or due to the fact that the answer is far form the global minimum. Local minimization techniques are also naturally sensitive to an initial estimate.

[0009] In general, a labeling f is a local minimum of the energy E if

[0010] E(f)≦E(f′) for any f′ “near to” f.  (1)

[0011] In case of discrete labeling, the labelings near to f are those that lie within a single move of f. Many local optimization techniques use standard moves, where only one pixel can change its label at a time. For standard moves, the equation above can be read as follows: if you are at a local minimum with respect to standard moves, then you cannot decrease the energy by changing a single pixel's label. In fact, this is a very weak condition. As a result, optimization schemes using standard moves frequently generate low quality solutions.

[0012] An example of a local method using standard moves is Iterated Conditional Modes (ICM). For each site (pixel or voxel), the label which gives the largest decrease of the energy function is chosen, until the iteration converges to a local minimum.

[0013] Another example of an algorithm using standard moves is simulated annealing. Simulated annealing randomly proposes some change in the state of the system. If the change results in a decrease of energy (which is equivalent to a decrease in cost in the more general sense of optimization), the change will always be taken. If it results in an increase in energy, it will be chosen using a probability scheme. At high temperatures (i.e. early in the simulated annealing process), changes of state that result in large increases in energy will be accepted with a higher probability than they would be at low temperatures (late in the simulated annealing process). Simulated annealing is widely used because it can optimize an arbitrary energy function. Minimizing an arbitrary energy function requires exponential time, and as a consequence simulated annealing is very slow.

[0014] Simulated annealing is inefficient partly because at each step, annealing changes the value of a single pixel. Theoretically, simulated annealing should eventually find the global minimum if it is run long enough. As a practical matter, it is necessary to decrease the algorithm's temperature parameter faster than required by the theoretically optimal schedule. Once annealing's temperature parameter is sufficiently low, the algorithm will converge to a local minimum with respect to standard moves.

[0015] An alternative method is to seek a local minimum using variational techniques, for example. Variational methods can be applied if the energy minimization problem is phrased in continuous terms. Variational techniques use Euler equations, which are guaranteed to hold at a local minimum. In continuous cases, the labels near to f in the equation above are normally defined as ∥f−f′∥≦ε where ε is a positive, constant and ∥·∥ is a norm, e.g. L₂over some appropriate functional space. To apply these algorithms to actual imagery requires discretization.

[0016] Another alternative is to use discrete relaxation labeling methods. In relaxation labeling, combinatorial optimization is converted into continuous optimization with linear constraints. Then some form of gradient descent, which gives the solution satisfying the constraints, is used.

[0017] There are also methods that have optimality guarantees in certain cases. Continuation methods, such as graduated non-convexity are an example. These methods involve approximating an intractable (non-convex) energy function by a sequence of energy functions beginning with a tractable (convex) approximation. There are circumstances where these methods are known to compute the optimal solution. Continuation methods can be applied to a large number of energy functions, but except for special cases, nothing is known about the quality of their output.

[0018] Mean field annealing is another popular minimization approach. It is based on estimating the partition function from which the minimum of the energy can be deduced. Computing the partition function, however, is computationally intractable, and saddle point approximations are used.

[0019] There are a few interesting energy functions where the global minimum can be rapidly computed via dynamic programming. Dynamic programming, however, is restricted essentially to energy functions in one-dimensional settings. In general, the two-dimensional energy functions that arise in early vision cannot be solved efficiently via dynamic programming.

[0020] It is an object of the present invention to provide a method and apparatus to improve the early vision stage of computer vision.

[0021] It is another object of the present invention to provide a method and apparatus to improve the speed while maintaining accuracy of minimizing energy functions.

[0022] The problems of early vision and minimizing energy functions are solved by the present invention of a method and apparatus for fast approximate energy minimization via graph cuts.

[0023] Many tasks in computer vision involve assigning a label, such as disparity for depth of field data to every pixel. Energy minimization may be used to accomplish this labeling. The present invention provides an efficient way of minimizing energy functions in assigning labels to pixels.

[0024] The major restriction is that the energy function's smoothness term must involve only pairs of pixels. Two methods of using graph cuts to compute a local minimum are described. They may be used even when very large moves are allowed.

[0025] The first move is an α-β swap. For a pair of labels, α, β, this move swaps the labels between an arbitrary set of pixels labeled a and another arbitrary set of pixels labeled β.

[0026] The first method generates a labeling such that there is no swap move that decreases the energy.

[0027] The second move is the α-expansion. For a label α, this move assigns an arbitrary set of pixels with the label α.

[0028] The second method which requires the smoothnes term to be a metric, generates a labeling such that there is no expansion move that decreases the energy.

[0029] The present invention together with the above and other advantages may best be understood from the following detailed description of the embodiments of the invention illustrated in the drawings, wherein:

[0030] FIG. 1ais a first algorithm for approximate minimization of energy functions as applied to early vision according to principles of the invention;

[0031] FIG. 1bis a second algorithm for approximate minimization of energy functions according to principles of the invention;

[0032] FIG. 2 is a graph system for the α-β swap method according to principles of the invention;

[0033] FIG. 3 are graph cut properties for the α-β swap method according to principles of the invention;

[0034] FIG. 4 is a flow chart of the operation of the α-β swap method of FIG. 1a;

[0035] FIG. 5 is a graph system for the α-expansion method of the present invention;

[0036] FIG. 6 are graph cut properties for the α-expansion method of FIG. 1b;

[0037] FIG. 7 is a flow chart of the operation of the α-expansion method of FIG. 1b;

[0038] FIGS. 8a, 8band 8care tables illustrating a local minimum where swap moves are arbitrarily far from the global minimum;

[0039] FIG. 9ais a graph system using the Potts model;

[0040] FIG. 9bis the graph system of FIG. 9aafter a multi-way cut method has been applied;

[0041] FIGS. 10a, 10b, 10c, 10d, 10e, and 10fshow the results from a real stereo pair with known ground truth;

[0042] FIG. 11 shows the graph of E_smoothversus time for the algorithms of the present invention versus simulated annealing;

[0043] FIG. 12 is a table giving the errors made by the expansion algorithm for different choices of λ;

[0044] FIG. 13ais the left image of a second stereo pair, FIG. 13bis the image of FIG. 13awhen minimizing E_P;

[0045] FIG. 13cis the image of FIG. 13awhen minimizing E_L;

[0046] FIG. 14ais a first image of a motion sequence where a cat moves against a moving background;

[0047] FIG. 14bshows the horizontal motions detected with the swap algorithm of the present invention in FIG. 14a;

[0048] FIG. 14cshows the vertical motions detected with the swap algorithm of the present invention in FIG. 14a;

[0049] FIG. 15ais another image; and

[0050] FIG. 15bshows the horizontal motions detected in the image of FIG. 15ausing the present invention.

[0051] FIGS. 1aand 1bare energy minimization functions of the present invention. Computer vision problems can be formulated in terms of energy minimization. The goal is to find solutions that are smooth over the surface of objects in a view but these solutions must also preserve sharp discontinuities at, for example, object boundaries.

[0052] In the algorithms of the present invention, an arbitrary initial solution is assumed. The initial solution may be completely arbitrary or it could be some less accurate solution to the computer vision problem then the present invention.

[0053] Energy functions are difficult and slow to solve. Approximate solutions provide good results in a reasonable period of time. The algorithms of the present invention provide approximate solutions with provable qualities. The algorithms of the present invention use graph cuts to obtain solutions. Each pixel in the view is a node. The graph can be thought of as a net and on the net are nodes and strings that connect the nodes. In a preferred embodiment of the invention, nodes are connected to neighboring nodes immediately above, below, to the right and to the left. The lines connecting the nodes represent the relationship between the nodes. In graph terminology, these lines are also referred to as “edges.” Each line belongs to a distinct pair of pixels. This node and line system is also referred to as the “neighborhood system.” In alternative embodiments of the invention, the graph could have diagonal links between nodes in addition to the links described above or instead of those links. Nodes could also be connected in ways other than to direct neighbors. Many types of neighborhood systems are possible and the invention is not limited to the neighborhood systems described herein.

[0054] V_p,qrepresents how pixels act within a pair of pixels. It is desirable for the difference between interacting pixels to be as small as possible. V_p,qcan also be thought of as a penalty for discontinuity between pixels. The present system puts restrictions on V_p,q. V_p,qmust be non-negative but it can be different for different pixel pairs. The penalties and restrictions are also part of the neighborhood system.

[0055] The input to the algorithms then is the neighborhood system, i.e. nodes, edges, penalty functions and assumptions.

[0056] The goal is to find a labeling f that minimizes the energy

[0057] E(f)=E_smooth(f)+E_data(f)  (2)

[0058] where E_smoothmeasures the extent to which f is not piecewise smooth, and

[0059] E_datameasures the disagreement between f and the observed data. The form of E_datameasures the disagreement between f and the observed data. The form of E_datais typically, $\begin{array}{cc}{E}_{\mathrm{data}}\left(f\right)=\sum _{p\in P}{D}_p\left(f_p\right)& \left(3\right)\end{array}$

[0060] where D_pmeasures how appropriate a label is for the pixel p given the observed data. In image restoration, for example, D_p(f_p) is typically (f_p−i_p), where i_pis the observed intensity of the pixel p.

[0061] The choice of E_smoothis important. There are many potential choices. For example, in regularization-based vision, E_smoothmakes f smooth everywhere. This gives poor results at object boundaries in general. Energy functions that do not have this problem are called discontinuity-preserving.

[0062] The major difficulty with energy minimization for early vision lies in the enormous computational costs. Typically these energy functions have many local minimums (i.e. they are non-convex). Also, the space of the possible labelings has dimension |P|, which is many thousands.

[0063] The energy functions used in the present invention are of the form, $\begin{array}{cc}E\left(f\right)=\sum _{\left\{p,q\right\}\in N}{V}_{p,q}\left(f_p,f_q\right)+\sum _{p\in P}{D}_p\left(f_p\right)& \left(4\right)\end{array}$

[0064] where N is the set of interacting pairs of pixels. Typically N consists of adjacent pixels, but it can be arbitrary. D_pis nonnegative but otherwise arbitrary. Only pairs of pixels interact directly in E_smooth. Each pair of pixels (p,q) can have its own distinct interaction penalty V_p,qindependent of any other pair of interacting pixels. To simplify notation, V_p,qwill be written as V.

[0065] The methods of the present invention approximately minimize the energy E(f) for an arbitrary finite set of labels L under two fairly general classes of interaction penalty V: metric and semi-metric: V is called a metric on the space of labels L if it satisfies:

[0066] V(α-β)=0⇄α=β  (5)

[0067] V(α-β)=V(α,β)≧0  (6)

[0068] V(α-β)≦V(α,γ)+V(β,γ)  (7)

[0069] for any labels α,β,γ∈L. If V does not satisfy the triangle inequality of equation (7), it is called semi-metric.

[0070] Both semi-metrics and metrics include important cases of discontinuity preserving interaction penalties. A penalty V is discontinuity preserving if it is bounded, i.e. sup_(x,y)∈RV(x,y)<∞. Examples of such interaction penalties for a one-dimensional label set L include the truncated quadratic V(α,β)=min(K,|α−β|²) (a semi-metric) and the truncated absolute distance V(α,β)=min(K,|α−β|) (a metric), where K is some constant. If L is multidimensional, |•| can be replaced by any norm. Another important discontinuity preserving interaction penalty is given by the Potts model V(α,β)=K·T(α≠β) (a metric), where T(•) is 1 if its argument is true, and otherwise 0.

[0071] The methods of the present invention generate a local minimum with respect to very large moves, i.e. more than one pixel at a time. The algorithms generate a labeling that is a local minimum of the energy for two types of large moves: α-expansion and α-β-swap. In contrast to the standard moves described above in the Background of the Invention, these moves allow large numbers of pixels to change their labels simultaneously. This makes the set of labelings within a single move of a locally optimal f exponentially large. For example, α-expansion moves are so strong that any labeling locally optimal with respect to these moves are within a known factor of the global minimum.

[0072] Any labeling f can be uniquely represented by a partition of image pixels P={P_l|l∈L} where P_l={p∈P|f_p=l} is a subset of pixels assigned label l. Since there is an obvious one to one correspondence between labelings f and partitions P, these notions will be used interchangingly.

[0073] Given a pair of labels α,β, a move from a partition P (labeling f) to a new partition P′ (labeling f′) is called an α-β swap if P_l=P′_lfor any label l≠α,β. This means that the only difference between P and P′ is that some pixels that were labeled α in P are labeled a in P′. A special case of an α-β swap is a move that gives the label α to some set of pixels that used to be labeled β.

[0074] Given a label α, a move from a partition P (labeling f) to a new partition P′ (labeling f′) is called an α-expansion if P_α⊂P′_αand P′_l⊂P_lfor any label l≠α. In other words, an α-expansion move allows any set of image pixels to change their labels to α.

[0075] A move which assigns a given label α to a single pixel is both an α-β swap and an α-expansion. As a consequence, a standard move is a special case of both an α-β swap and an α-expansion.

[0076] Algorithms and Properties

[0077] Referring again to FIG. 1, there are two minimization algorithms in the present invention. The swap algorithm of FIG. 1afinds a local minimum when swap moves are allowed, and the expansion algorithm of FIG. 1bfinds a local minimum when the expansion moves are allowed. Finding such a local minimum is not a trivial task. Given a labeling f, there is an exponential number of swap and expansion moves. Therefore, even checking for a local minimum requires exponential time if it is performed naively. In contrast checking for a local minimum when only the standard moves are allowed is easy since there is only a linear number of standard moves given any labeling f.

[0078] The algorithms are efficient graph based methods to find the optimal α-β or α-expansion given a labeling f. This is the key step in the algorithms. Once these methods are available, it is easy to design some variations of the “fastest descent” technique that can efficiently find the corresponding local minimums.

[0079] The structure of the algorithms is similar. A single execution of steps 3.1-3.2 is called an iteration, and an execution of steps 2-4 is a cycle. In each cycle, the algorithm performs an iteration for every label (expansion algorithm) or for every pair of labels (swap algorithm), in a certain order that can be fixed or random. A cycle is successful if a strictly better labeling is found at any iteration. The algorithms stop after the first unsuccessful cycle since no further improvement is possible. Obviously, a cycle in the swap algorithm takes |L|²iterations, and a cycle in the expansion algorithm takes |L| iterations. These algorithms are guaranteed to terminate in a finite number of cycles. Termination can be proved to be in O(|P|) cycles. In experiments, the algorithm stops after a few cycles, and most of the improvements occur during the first cycle.

[0080] Graph cuts are used to efficiently find {circumflex over (f)} for the key part of each algorithm in step 3.1. Step 3.1 uses a single minimum cut on a graph whose size is O(|P|). The graph is dynamically updated after each iteration.

[0081] Graph Cuts

[0082] Before describing the key step 3.1 of the swap and the expansion algorithms, graph cuts will be reviewed. Let G=(V,ε) be a weighted graph with two distinguished vertices called the terminals. A cut C⊂ε is a set of edges such that the terminals are separated in the induced graph G(C)=(V,ε−C). In addition, no proper subset of C separates the terminals in G(C). The cost of the cut C, denoted |C|, equals the sum of its edge weights.

[0083] The minimum cut problem is to find the cut with smallest cost. This problem can be solved efficiently by computing the maximum flow between the terminals, according to a theorem due to Ford and Fulkerson. There are a large number of fast algorithms available in the art for this problem. The worst case complexity is low-order polynomial; however, for graphs of the present invention, the running time is nearly linear in practice.

[0084] Finding the Optimal Swap Move

[0085] Given an input labeling f (partition P) and a pair of labels α, β, a labeling {circumflex over (f)} needs to be found that minimizes E over all labelings within one α-β swap of f. This is the, critical step in the swap move algorithm of FIG. 1a. The other technique of the present invention is based on computing a labeling corresponding to a minimum cut on a graph G_αβ=(V_αβ,ε_αβ). The structure of this graph is dynamically determined by the current partition P and by the labels α,β.

[0086] This section is organized as follows. The construction of G_αβ for a given f (or P) is described. Cuts C on G_αβ correspond in a natural way to labelings f^Cwhich are within one α-β swap move of f. Theorem 4.4 shows that the cost of a cut is |C|=E(f^C) plus a constant. A corollary from this theorem states the main result that the desired labeling {circumflex over (f)} equals f^Cwhere C is a minimum cut on G_αβ.

[0087] FIG. 2 shows the structure of the graph. For legibility, this figure shows the case of a 1-D image. For any image, the structure of G_αβ100 will be as follows. The set of vertices includes the two terminals α 102 and β 104, as well as image pixels P 106 in the sets P_α and P_β (that is f_p∈{α,β}). “α” and “β” are also referred to as labels. Thus, the set of vertices V_αβ consists of α, β, and P_αβ=P_α∪P_β. Each pixel p∈P_αβ is connected to the terminals α and β by edges t_p^α108 and t_p^β110, respectively. For brevity, these edges are referred to as t-links (terminal links). Each pair of pixels {p,q}⊂P_αβ which are neighbors (i.e. {p,q}∈N is connected by an edge e_{p,q}112 which are called an n-link (neighbor link). The set of edges ε_αβ thus consists of U_p∈P{t_p^α,t_p^β} (the t-links) and $U_{\begin{array}{c}\left\{p,q\right\}\in N\\ p,q\in P_{\mathrm{\alpha \beta }}\end{array}}{e}_{\left\{p,q\right\}}$

[0088] (the n-links). The weights assigned to the edges according to Table 1.

[0089] Any cut C on G_αβ must sever (include) exactly one t-link for any pixel p∈P_αβ. If neither t-link were in C, there would be a path between the terminals; while if both t-links were cut, then a proper subset of C would be a cut. Thus, any cut leaves each pixel in P_αβ with exactly one t-link. This defines a natural labeling f^Ccorresponding to a cut C on G_αβ, $\begin{array}{cc}{f}_p^C\{\begin{array}{cc}\alpha & \mathrm{if}{t}_p^{\alpha }\in C\mathrm{for}p\in P_{\mathrm{\alpha \beta }}\\ \beta & \mathrm{if}{t}_p^{\beta }\in C\mathrm{for}p\in P_{\mathrm{\alpha \beta }}\\ f_p& \mathrm{for}p\in P,p\notin P_{\mathrm{\alpha \beta }}\end{array}& \left(8\right)\end{array}$

[0090] In other words, if the pixel p is in P_αβ hen p is assigned label α when the cut C separates p from the terminal α; similarly, p is assigned label β when C separates p from the terminal β. If p is not in P_αβ then its initial label f_pis retained. This implies

[0091] Lemma 4.1 A labeling f^Ccorresponding to a cut C on G_αβ is one β-β swap away from the initial labeling f.

[0092] It is easy to show that a cut C severs an n-link e_{p,q} between neighboring pixels on G_αβ if and only if C leaves the pixels p and q connected to different terminals. FIG. 3 shows Property 4.2, which for any cut C and for any n-link e_{p,q} are:

[0093] a) If t_p^α,t_p^α∈Cthen e_{p,q}∉C.

[0094] b) If t_p^β,t_q^β∈Cthen e_{p,q}∉C.

[0095] c) If t_p^β,t_q^α∈Cthen e_{p,q}∈C.

[0096] d) If t_p^α,t_q^β∈Cthen e_{p,q}∈C.  (9)

[0097] Properties (a) and (b) follow from the requirement that no proper subset of C should separate the terminals. Properties (c) and (d) also use the fact that a cut has to separate the terminals.

[0098] The next lemma is a consequence of property 4.2 and equation 8.

[0099] Lemma 4.3 For any cut C and for any n-link e_{p,q}

[0100] |C∩e_{p,q}|=V(f_p^C,f_q^C)

[0101] PROOF: There are four cases with similar structure; the case where t_β^α,t_q^β∈C will be illustrated. In this case, e_{p,q}∈C and, therefore, |C∩ε_{p,q}|=V(α,β). By (8), f_p^C=α and f_q^eα=β.

[0102] Note that this proof assumes that V is a semi-metric, i.e. that equations 5 and 6 hold.

[0103] Lemmas 4.1 and 4.3 plus property 4.2 yield

[0104] Theorem 4.4 There is a one to one correspondence between cuts C on G_αβ and labelings hat are one α-β swap from f. Moreover, the cost of a cut C on G_αβ is |C|=E(f^C) plus a constant.

[0105] Proof: The first part follows from the fact that the severed t-links uniquely determine the labels assigned to pixels p and the n-links that must be cut. The cost of a cut C is $\begin{array}{cc}\left|C\right|=\sum _{p\in P_{\mathrm{\alpha \beta }}}^{}\left|C\bigcap \left\{t_p^{\alpha },t_p^{\beta }\right\}\right|+\sum _{\begin{array}{c}\left\{p,q\right\}\in N\\ \left\{p,q\right\}⋐P_{\mathrm{\alpha \beta }}\end{array}}^{}\left|C\bigcap e_{\left\{p,q\right\}}\right|& \left(10\right)\end{array}$

[0106] Note that for p∈P_αβ, $\begin{array}{cc}\left|C\bigcap \left\{t_p^{\alpha },t_p^{\beta }\right\}\right|=\{\begin{array}{c}\left|t_p^{\alpha }\right|\mathrm{if}{t}_p^{\alpha }\in C\\ \left|t_p^{\beta }\right|\mathrm{if}{t}_p^{\alpha }\in C\end{array}=D_p\left(f_p^c\right)+\sum _{\begin{array}{c}q\in N_p\\ q\notin P_{\mathrm{\alpha \beta }}\end{array}}^{}V\left(f_p^C,f_q\right)& \left(11\right)\end{array}$

[0107] Lemma 4.3 gives the second term in (10). Thus, the total cost of a cut C is $\begin{array}{cc}\begin{array}{c}\left|C\right|=\sum _{p\in P_{\mathrm{\alpha \beta }}}^{}{D}_p\left(f_p^C\right)+\sum _{p\in P_{\mathrm{\alpha \beta }}}^{}\sum _{\begin{array}{c}q\in N_p\\ q\notin P_{\mathrm{\alpha \beta }}\end{array}}^{}V\left(f_p^c,f_q\right)+\sum _{\begin{array}{c}\left\{p,q\right\}\in N\\ \left\{p,q\right\}⋐P_{\mathrm{\alpha \beta }}\end{array}}^{}V\left(f_p^c,f_q\right)\\ =\sum _{p\in P_{\mathrm{\alpha \beta }}}^{}{D}_p\left(f_p^C\right)+\sum _{\begin{array}{c}\left\{p,q\right\}\in N\\ p\mathrm{or}q\in P_{\mathrm{\alpha \beta }}\end{array}}^{}V\left(f_p^C,f_q^C\right)\end{array}& \left(12\right)\end{array}$

[0108] This can be rewritten as |C|=E(f^C)−K where $\begin{array}{cc}K=\sum _{p\notin P_{\mathrm{\alpha \beta }}}^{}{D}_p\left(f_p\right)+\sum _{\begin{array}{c}\left\{p,q\right\}\in N\\ \left\{p,q\right\}\bigcap P_{\mathrm{\alpha \beta }=0}\end{array}}^{}V\left(f_p^{},f_q^{}\right)& \left(13\right)\end{array}$

[0109] is the same constant for all cuts C.

[0110] Corollary 4.5 The lowest energy labeling within a single α-β swap move from f is {circumflex over (f)}=f^C, where C is the minimum cut on G_αβ.

[0111] FIG. 4 shows the operation of the application of the α-β swap method. A single view typically has a plurality of labels. The input to the method, as described above, is the neighborhood system with some arbitrary or other initial solution having an initial labeling, block 150.

[0112] Next, a pair of labels to swap is chosen, block 152. Then, using the graph cuts described above, the best swap for the chosen pair of labels is determined, block 154. The minimum cut is chosen by assigning weights to the various links in the graph according to Table 1. A minimum cut is calculated from the weighted links.

[0113] The new labeling that results from finding the best swap move is stored, block 156.

[0114] In the preferred embodiment of the invention, the method determines if any pair of label combinations can be improved, block 158. If yes, the algorithm moves on to another pair of labels and the process repeats until no label combinations can be improved. If there are no more improvements that can be made, the process ends, block 160. In alternative embodiments of the invention, the iteration continues only until some acceptable solution is obtained.

[0115] Finding the Optimal Expansion Move

[0116] Given an input labeling f (partition P) and a label α, a labeling {circumflex over (f)} that minimizes E over all labelings within one α-expansion of f needs to be found. This is the critical step in the expansion move algorithm of FIG. 1b. A technique that solves the problem assumes that each V is a metric, and thus satisfies the triangle inequality (7). The technique is based on computing a labeling corresponding to a minimum cut on a graph G_α=(V_α,ε_α). The structure of this graph is determined by the current partition P and by the label α. As before, the graph dynamically changes after each iteration.

[0117] This section is organized as follows. The construction of G_α for a given f (or P) and a is described. Then cuts C on G_α are shown to correspond in a natural way to labelings f^Cwhich are within one α-expansion move of f. Then, based on a number of simple properties, a class of elementary cuts is defined. Theorem 5.4 shows that elementary cuts are in one to one correspondence with labelings that are within one α-expansion of f, and also that the cost of an elementary cut is |C|=E(f^C). A corollary from this theorem states he main result that the desired labeling {circumflex over (f)} is f^Cwhere C is a minimum cut on G_α.

[0118] The structure of the graph is illustrated in FIG. 5. For legibility, this figure shows the case of a 1-D image. The set of vertices includes the two terminals a 202 and {overscore (α)} 204, as well as all image pixels p∈P 206. In addition, for each pair of neighboring pixels {p,q}∈N separated in the current partition (i.e. such that f_p≠f_q, an auxiliary vertex a_{p,q} is created. Auxiliary nodes are introduced at the boundaries between partition sets P_lfor l∈L. Thus, the set of vertices is $\begin{array}{cc}{V}_{\alpha }=\left\{\alpha ,\stackrel{\_}{\alpha },P,\bigcup _{\begin{array}{c}\left\{p,q\right\}\in N\\ f_p\ne f_q\end{array}}{a}_{\left\{p,q\right\}}\right\}& \left(14\right)\end{array}$

[0119] Each pixel p∈P is connected to the terminals α and {overscore (α)} by t-links t_p^α208 and t_p^{{overscore (α)}}210, respectively. Each pair of neighboring pixels {p,q}∈N which are not separated by the partition P (i.e. such that f_p=f_q) is connected by an n-link; e_{p,q}212. For each pair of neighboring pixels {p,q}∈N such that f_p≠f_qa triplet of edges ε_{p,q}={e_{p,a},e_{a,q},t_a^{{overscore (α)}}} is created where a=a_{p,q} is the corresponding auxiliary node. The edges e_{p,a} and e_{a,q} connect pixels p and q to a_{p,q} and the t-link t_a^{{overscore (α)}} connects the auxiliary node a_{p,q} to the terminal {overscore (α)}. The set of all edges may be written as $\begin{array}{cc}{ϵ}_{\alpha }=\left\{\bigcup _{p\in P}\left\{t_p^{\alpha },t_p^{\stackrel{\_}{\alpha }}\right\},\bigcup _{\begin{array}{c}\left\{p,q\right\}\in N\\ f_p\ne f_q\end{array}}{ϵ}_{\left\{p,q\right\}},\bigcup _{\begin{array}{c}\left\{p,q\right\}\in N\\ f_p=f_q\end{array}}{e}_{\left\{p,q\right\}}\right\}& \left(15\right)\end{array}$

[0120] The weights are assigned to the edges according to, Table 2.

[0121] Any cut C on G_αmust sever (include) exactly one t-link for any pixel p∈P. This defines a natural labeling f^Ccorresponding to a cut C on G_α. Formally, $\begin{array}{cc}{f}_p^c=\{\begin{array}{cc}\alpha \mathrm{if}{t}_p^{\alpha }\in C& \\ & {\forall }_p\in P\\ f_p\mathrm{if}{t}_p^{\stackrel{\_}{\alpha }}\in C& \end{array}& \left(16\right)\end{array}$

[0122] In other words, a pixel p is assigned label α if the cut C separates p from the terminal α, while p is assigned its old label f_pif C separates p from {overscore (α)}. Note that for p∉P_α the terminal {overscore (α)} represents labels assigned to pixels in the initial labeling f. Clearly this leads to:

[0123] Lemma 5.1 A labeling f^Ccorresponding to a cut C on G_α is one α-expansion away from the initial labeling f.

[0124] It is also easy to show that a cut C severs an n-link e_{p,q} between neighboring pixels {p,q}∈N such that f_p=f_qif and only if C leaves the pixels p and q connected to different terminals. In other words, property 4.2 holds when “{overscore (α)}” is substituted for “β”. This as property is referred to as 4.2({overscore (α)}). Analogously, that property 4.2 and equation 16 establish lemma 4.3 for the n-links e_{p,q} in G_α.

[0125] Consider now the set of edges ε_{p,q} corresponding to a pair of neighboring pixels {p,q}∈N such that f_p≠f_q. In this case, there are several different ways to cut these edges even when the pair of severed t-links at p and q is fixed. However, a minimum cut C on G_α is guaranteed to sever the edges in ε_{p,q} depending on what t-links are cut at the pixels p and q.

[0126] FIG. 6 shows the rule for this case as described in property 5.2 below. Assume that α=α_{p,q} is an auxiliary node between the corresponding pair of neighboring pixels.

[0127] Property 5.2 If {p,q}∈N and f_p≠f_qthen a minimum cut C on G_α satisfies:

[0128] a) If t_p^α,t_q^α∈Cthen C∩ε_{p,q}=0.

[0129] b) If t_p^{{overscore (α)}},t_q^{{overscore (α)}}∈Cthen C∩ε_{p,q}=t_a^{{overscore (α)}}.

[0130] c) If t_p^{{overscore (α)}},t_q^α∈Cthen C∩ε_{p,q}=e_{p,a}

[0131] d) If t_p^α,t_q^{{overscore (α)}}∈Cthen C∩ε_{p,q}=e_{a,q}

[0132] Property 5.2(a) results from the fact that no subset of C is a cut. The others follow from the minimality of |C| and the fact that |e_{p,a}|, |e_{a,q}| and |t_a^{{overscore (α)}} satisfy the triangle inequality so that cutting any one of them is cheaper than cutting the other two together.

[0133] Lemma 5.3 If If {p,q}∈N and f_p≠f_qthen the minimum cut C on G_α satisfies:

[0134] |C∩ε_{p,q}|=V(f_p^C,f_q^C)

[0135] Proof: The equation follows from property 5.2, equation (16), and the edge weights. For example,

[0136] If t_p^{{overscore (α)}},t_q^{{overscore (α)}}∈C then C∩ε_{p,q}=|t_a^{{overscore (α)}}=V(f_p,f_q) At the same time, (16) implies that f_p^C=f_pand f_p^C=f_q. Note that the right penalty V is imposed whenever f_p^C≠f_q^C, due to the auxiliary pixel construction.

[0137] Property 4.2 ({overscore (α)}) holds for any cut, and property 5.2 holds for a minimum cut. However, there can be other cuts; besides the minimum cut that satisfy both properties. An elementary cut on G_α is defined as a cut that satisfies properties 4.2 ({overscore (α)}) and 5.2.

[0138] Theorem 5.4 Let G_α be constructed as above given f and α. Then there is a one to one correspondence between elementary cuts on G_α and labelings within one α-expansion of f. Moreover, for any elementary cut C, |C|=E(f^C) is the result. Proof: An elementary cut C is uniquely determined by the corresponding labeling f^C. The label f_p^Cat the pixel p determines which of the t-links to p is in C. Property 4.2 ({overscore (α)}) shows which n-links e_{p,q} between pairs of neighboring pixels {p,q} such that f_p=f_qshould be severed. Similarly, property 5.2 determines which of the links in ε_{p,q} corresponding o {p,q}∈N such that f_p≠f_qshould be cut.

[0139] The cost of an elementary cut C is $\begin{array}{cc}\left|C\right|=\sum _{p\epsilon P}\left|C\bigcap \left\{t_p^{\alpha },t_p^{\stackrel{\_}{\alpha }}\right\}\right|+\sum _{\underset{f_p=f_q}{\left(p,q\right)}\epsilon N}\left|C\bigcap e_{\left\{p,q\right\}}\right|+\sum _{\underset{f_p\ne f_q}{\left(p,q\right)}\epsilon N}\left|C\bigcap {ϵ}_{\left\{p,q\right\}}\right|& \left(17\right)\end{array}$

[0140] It is easy to show that for any pixel p∈P, the result is

[0141] |C∩{t_p^α,t_p^{{overscore (α)}}}|=D_p(f_p^C).  (18)

[0142] Lemmas 4.3 and 5.3 hold for elementary cuts, since they were based on properties 4.2 and 5.2. Thus, the total cost of a elementary cut C is $\begin{array}{cc}\left|C\right|=\sum _{p\epsilon P}{D}_p\left(f_p^C\right)+\sum _{\left(p,q\right)\epsilon N}V\left(f_p^C,f_q^C\right)=E\left(f^C\right)& \left(19\right)\end{array}$

[0143] Therefore, |C|=E(f^C).

[0144] Our main result is a simple consequence of this theorem, since the minimum cut is an elementary cut.

[0145] Corollary 5.5 The lowest energy labeling within a single α-expansion move from f is {circumflex over (f)}=f^C, where C is the minimum cut on G_α.

[0146] FIG. 7 shows the operation of the α-expansion method of the present invention. A single view in computer vision typically has a plurality of labels. The input to the method is the neighborhood system with some arbitrary or other initial solution including a set of initial labeling, block 250.

[0147] Next, a label is chosen, block 252. The chosen label is expanded, block 254. That is, any pixel is allowed to take on the chosen label. Then the best solution is determined using the minimum graph cut method, block 256. Each link in the graph is assigned a weight according to Table 2 and cuts are made to determine the best solution as described above. In a preferred embodiment of the invention, the method then determines if any more improvements may be made, block 258. If yes, the method repeats until no further improvements may be made and then the process ends, block 260. In an alternative embodiment of the invention, the iteration continues only until an acceptable solution is obtained.

[0148] Optimality Properties

[0149] Any local minimum generated by the expansion moves algorithm of the present invention is within a known factor from the global optimum. This algorithm works in case of metric V. The swap move algorithm can be applied to a wider class of semi-metric V's but, unfortunately, it does not have any guaranteed optimality properties. A provably good solution can be obtained even for semi-metric V by approximating such V's with a simple Potts metric.

[0150] The Expansion Move Algorithm

[0151] A local minimum when expansion moves are allowed is within a known factor of the global minimum. This factor, which can be as small as 2, will depend on V. Specifically, let $\begin{array}{cc}c=\frac{{\max }_{\alpha \ne \beta \epsilon L}V\left(\alpha ,\beta \right)}{{\min }_{\alpha \ne \beta \epsilon L}V\left(\alpha ,\beta \right)}& \left(20\right)\end{array}$

[0152] be the ratio of the largest non zero value of V to the smallest non zero value of V. Note that c is well defined since V(α,β)≠0 for α≠β according to the metric properties (2) and (3). If V_p,q's are different for neighboring pairs p,q $\mathrm{then}c={\max }_{p,q\epsilon N}\left(\frac{{\max }_{\alpha \ne \beta \epsilon L}V\left(\alpha ,\beta \right)}{{\min }_{\alpha \ne \beta \epsilon L}V\left(\alpha ,\beta \right)}\right).$

[0153] Theorem 6.1 Let {circumflex over (f)} be a local minimum when the expansion moves are allowed and f* be the globally optimal solution. Then E({circumflex over (f)})≦2cE(f*)

[0154] Proof: Fix some α∈L and let

[0155] P_α={p∈P|f*_p=α}  (21)

[0156] A labeling f^α can be produced within one α-expansion move from {circumflex over (f)} as follows: $\begin{array}{cc}{f}_p^{\alpha }=\{\begin{array}{cc}\alpha & \mathrm{if}p\in P_{\alpha }\\ {\hat{f}}_p& \mathrm{otherwise}\end{array}& \left(22\right)\end{array}$

[0157] The key observation is that since {circumflex over (f)} is a local minimum if expansion moves are allowed,

[0158] E({circumflex over (f)})≦E(f^α).  (23)

[0159] Let S be a set consisting of any number of pixels in P and any number of pairs of neighboring pixels in N. E(f|S) is defined to be a restriction of the energy of labeling f to the set S: $\begin{array}{cc}E(f\left|S\right)=\sum _{p\epsilon S}{D}_p\left(f_p\right)+\sum _{\left(p,q\right)\epsilon S}V\left(f_p,f_q\right)& \left(24\right)\end{array}$

[0160] Let I^α be the set of pixels and pairs of neighboring pixels contained inside P_α. Also, let B^α be the set of pairs of neighboring pixels on the boundary of P_α and O^α be the set of pixels and pairs of neighboring pixels contained outside of P_α. Formally, $\begin{array}{cc}{I}^{\alpha }=P_{\alpha }\bigcup \left\{\left\{p,q\right\}\in N:p\in P_{\alpha },q\in P_{\alpha }\right\}& \left(25\right)\\ B^{\alpha }=\left\{\left\{p,q\right\}\in N:p\in P_{\alpha },q\notin P_{\alpha }\right\}& \left(26\right)\\ O^{\alpha }=\left(P-P_{\alpha }\right)\bigcup \left\{\left\{p,q\right\}\in N:p\notin P_{\alpha },q\notin P_{\alpha }\right\}& \left(27\right)\end{array}$

[0161] The following three facts hold:

[0162] E(f^α|O^α)=E({circumflex over (f)}|O^α),  (28)

[0163] E(f^α|I^α)=E(f*|I^α),  (29)

[0164] E(f^α|B^α)≦E(f*|B^α).  (30)

[0165] Equations 28 and 29 are obvious from the definitions in (22) and (21). Equation (30) holds because for any {p,q}∈B^α there is V(f_p^α,f_q^α)≦cV(f*_p,f*_q)≠0.

[0166] Since I^α∪B^α∪O^α includes all pixels in P and all neighboring pairs of pixels in N, both sides of (23) can be expanded to get:

[0167] E({circumflex over (f)}|I^α)+E({circumflex over (f)}|B^α)+E({circumflex over (f)}|O^α)≦E(f^α|I^α)+E(f^α|B^α)+E(f^α|O^α)  (31)

[0168] Using (28), (29) and (30) the following equation is derived from the equation above:

[0169] E({circumflex over (f)}|I^α)+E({circumflex over (f)}|B^α)≦E(f*|I^α)+cE(f*|B^α)  (32)

[0170] To get the bound on the total energy, equation (32) is summed over all labels α∈L: $\begin{array}{cc}\sum _{\alpha \epsilon L}(E(\hat{f}\left|I^{\alpha }\right)+E\left(\hat{f}\left|B^{\alpha }\right)\right)\le \sum _{\alpha \epsilon L}\left(E\left(f^{*}❘I^{\alpha }\right)+\mathrm{cE}\left(f^{*}❘B^{\alpha }\right)\right)& \left(33\right)\end{array}$

[0171] Let $B=\bigcup {\underset{\alpha \epsilon L}{B}}^{\alpha }.$

[0172] Observe that for every {p,q}∈B, the term V({circumflex over (f)}_p,{circumflex over (f)}_q)=E({circumflex over (f)}|{p,q}) appears twice on the left side of (33), once in E({circumflex over (f)}|B^α) for α=f*_pand once in E({circumflex over (f)}|B^α) for α=f*_q. Similarly every V(f*_p,f*_q*)=E(f*|{p,q}) appears 2c times on the right side of (33). Therefore equation (33) can be rewritten to get the bound of 2c:

[0173] E({circumflex over (f)})+E({circumflex over (f)}|B)≦E(f*)+(2c−1)E_B(f*)≦2cE(f*)  (34)

[0174] It can be shown that any algorithm that is within a factor of 2 for the Potts model is within a factor of 2c for an arbitrary metric V.

[0175] Approximating a Semi-metric

[0176] FIGS. 8a, 8band 8care tables illustrating a local minimum where swap moves are arbitrarily far from the global minimum.

[0177] The expansion algorithm of the present invention can be used to get an answer within a factor of 2c from the optimum of energy (4) even when V is a semi-metric. Actually just V(α,β)≧0 and V(α,β)=0⇄α=β are needed. Here c is the same as in Theorem 6.1. This c is still well defined for a semi-metric. Suppose that penalty V inside the definition of energy E in (4) is a semi-metric. Let r be any real number in the interval [m,M] where $m=\underset{\alpha \ne \mathrm{\beta \epsilon }L}{\min }V\left(\alpha ,\beta \right)\mathrm{and}M=\underset{\alpha \ne \mathrm{\beta \epsilon }L}{\max }V\left(\alpha ,\beta \right)$

[0178] Define a new energy based on the Potts interaction model $\begin{array}{cc}{E}_P\left(f\right)=\sum _{p\epsilon P}{D}_p\left(f_p\right)+\sum _{\left\{p,q\right\}\epsilon N}r·T\left(f_p\ne f_q\right)& \left(35\right)\end{array}$

[0179] Theorem 6.2 If {circumflex over (f)} is a local minimum of E_Pgiven the expansion moves and f* is the global minimum of E(f) then E({circumflex over (f)})≦2cE(f*).

[0180] Proof: Suppose f° is the global minimum of E_P. Then $\frac{r}{M}E\left(\hat{f}\right)\le E_P\left(\hat{f}\right)\le 2E_p\left(f^o\right)\le 2E_p\left(f^{*}\right)\le 2\frac{r}{m}E\left(f^{*}\right)$

[0181] where the second inequality follows from Theorem 6.1. Note that c=M/m.

[0182] Thus to find an answer within a fixed factor from the global minimum for a semi-metric V, one can take a local minimum {circumflex over (f)} given the expansion moves for E_Pas defined above. Note that that such an {circumflex over (f)} is not a local minimum of E(f) given the expansion moves. In practice, it has been found find that local minimum given the swap moves gives empirically better results than using {circumflex over (f)}. In fact, the estimate {circumflex over (f)} can be used as a good starting point for the swap algorithm. In this case, the swap move algorithm will also generate a local minimum whose energy is within a known factor from the global minimum.

[0183] The Potts Model

[0184] An interesting special case of the energy in equation (4) arises when V is given by the Potts model $\begin{array}{cc}{E}_P\left(f\right)=\sum _{\left\{p,q\right\}\epsilon N}{u}_{\left\{p,q\right\}}·T\left(f_p\ne f_q\right)+\sum _{p\epsilon P}{D}_p\left(f_p\right)& \left(36\right)\end{array}$

[0185] In this case, discontinuities between any pair of labels are penalized equally. This is in some sense the simplest discontinuity preserving model and it is especially useful when the labels are unordered or the number of labels is small. The Potts interaction penalty V_p,q=u_{p,q}·T(f_p≠f_q) is a metric; in this case c=1 and the expansion algorithm of the present invention gives a solution that is within a factor of 2 of the global minimum. Note that by definition c≧1, so this is the energy function with the best bound.

[0186] The Potts model energy minimization problem is closely related to a known combinatorial optimization problem called the multiway cut problem. The Potts model energy minimization problem can be reduced to the multiway cut problem. More precisely, it can be proved that the global minimum of the Potts model energy E_Pcan be computed by finding the minimum cost multiway cut on an appropriately constructed graph. In order to efficiently compute the global minimum of E_P, a certain class of multiway cut problems that are known to be NP-hard could also be solved. Minimizing E_Pis NP-hard, and so is minimizing the energy in (4).

[0187] The multiway cut problem is defined on a graph G=(V,E), with non-negative edge weights, with a set of terminal vertices L⊂V. A subset of the edges C⊂ε is called a multiway cut if the terminals are completely separated in the {induced} graph G(C)=(V, ε−C). It is also required that no proper subset of C separates the terminals in G(C). The cost of the multiway cut C is denoted by |C| and equals the sum of its edge weights. The multiway cut problem is to find the minimum cost multiway cut. The multiway cut problem is NP-: complete. The multiway cut problem is a generalization of the standard two-terminal graph cut problem.

[0188] The Potts Model and the Multiway Cut Problem

[0189] The problem of minimizing the Potts energy E_P(f) can be solved by computing a minimum cost multiway cut on a certain graph. Given V=P∪L, means that G contains two types of vertices: p-vertices (pixels) and l-vertices (labels). Note that l-vertices will serve as terminals or the multiway cut problem. Two p-vertices are connected by an edge if and only if the corresponding pixels are neighbors in the neighborhood system N. The set ε_Nconsists of the edges between p-vertices, called n-links. Each n-link {p,q}∈ε_Nis assigned a weight w_{p,q}=u_{p,q}.

[0190] Each p-vertex is connected by an edge to each l-vertex. An edge {p,l} that connects a p-vertex with a terminal (an l-vertex) will be called a t-link and the set of all such edges will be denoted by ε_τ. Each t-link {p,l} in ε_τ is assigned a weight w_{p,l}=K_p−D_p(l); where K_p>max_lD_p(l) is a constant that is large enough to make the weights positive. The edges of the graph are ε=ε_N∪ε_τ.

[0191] FIG. 9ashows the structure of the graph G. The graph G=(V,ε) with terminals L={1, . . . , k}. The pixels p∈P are shown as white squares. Each pixel has an n-link to its four neighbors. Each pixel is also connected to all terminals by t-links (some of the t-links are omitted from the drawing for clarity). The set of vertices V=P∪L includes all pixels and terminals. The set of edges ε=ε_N∪ε_τ, consists of all n-links and t-links.

[0192] It is easy to see that there is a one-to-one correspondence between multiway cuts and labelings. A multiway cut C corresponds to the labeling f^Cwhich assigns the label l to all pixels p which are t-linked to the l-vertex in G(C). An example of a multiway cut and the corresponding image partition (labeling) is given in FIG. 9b. The multiway cut corresponds to a unique partition (labeling) of image pixels.

[0193] Theorem 7.1 If C is a multiway cut on G_α then |C|=E_P(f^C) plus a constant.

[0194] Corollary 7.2 If C is a minimum cost multiway cut on G, then f^Cminimizes E_P.

[0195] While the multiway cut problem is known to be NP-complete if there are more than 2 terminals, there is a fast approximation algorithm. This algorithm works as follows. First, for each terminal l∈L, it finds an isolating two-way minimum cut C(l) that separates l from all other terminals. This is just the standard graph cut problem. Then the algorithm generates a multiway cut C=U_l≠lC(l) where l_max=arg max_l∈L|C(l)| is the terminal with the largest cost isolating cut. This “isolation heuristic” algorithm produces a cut which is optimal to within a factor of $2-\frac{2}{\left|L\right|}.$

[0196] However, the isolation heuristic algorithm suffers from two problems that limits its applicability to the energy minimization problem.

[0197] First, the algorithm will assign many pixels a label that is chosen essentially arbitrarily. Note that the union of all isolating cuts U_l∈LC(l) may leave some vertices disconnected from any terminal. The multiway cut C=U_l≠lC(l) connects all those vertices to the terminal l_max.

[0198] Second, while the multiway cut C produced is close to optimal, this does not imply that the resulting labeling f^Cis close to optimal. Formally, write theorem 7.1 as |C|=E_P(C)+K (the constant K results from the K_p's. The isolation heuristic gives a solution Ĉ such that |Ĉ|≦2|C*|, where C* is the minimum cost multiway cut. Thus, E_P(Ĉ)+K≦2(E_P(C*)+K), so E_P(Ĉ)≦2E_P(C*)+K. As a result, the isolation heuristic algorithm does not produce a labeling whose energy is within a constant factor of optimal.

[0199] Minimizing the Potts Energy is NP-hard

[0200] For an arbitrary fixed graph G=(V,E), it is possible to construct an instance of minimizing E_P(f) where the optimal labeling f* determines a minimum multiway cut on G. This will prove that a polynomial-time method for finding f* would provide a polynomial-time algorithm for finding the minimum cost multiway cut, which is known to be NP-hard.

[0201] The energy minimization problem takes as input a set of pixels P, a neighborhood relation N and a label set L, as well as a set of weights u_{p,q} and a function D_p(l). The problem is to find the labeling f* that minimizes the energy E_P(f) given in equation (36).

[0202] Let G=(V,ε) be an arbitrary weighted graph with terminal vertices {t₁, . . . , t_k}⊂V and edge weights w_{p,q}. The energy minimization using P=V, N=E, and u_{p,q}=w_{p,q}. The label set will be L={1, . . . , k}. Let K be a constant such that K>E_P(f*); for example, K can be selected to be the sum of all w_{p,q}. The function D_p(l) will force f*(t_j)=j; if p=t_jis a terminal vertex, $D_P\left(l\right)=\{\begin{array}{cc}0& l=j,\\ K& \mathrm{otherwise}.\end{array}$

[0203] For a non-terminal vertex p all labels are equally good,

[0204] ∀l D_p(l)=0.

[0205] A labeling f is defined to be feasible if the set of pixels labeled j by f forms a connected component that includes t_j. Feasible labelings obviously correspond one-to-one with multiway cuts.

[0206] Theorem 7.3 The labeling f* is feasible, and the cost of a feasible labeling is the cost of the corresponding multiway cut.

[0207] Proof: To prove that f* is feasible, suppose that there were a set S of pixels that f* labeled j which were not part of the component containing t_j. A labeling could then be obtained having lower energy by switching this set to the label of some pixel on the boundary of S. The energy of a feasible labeling f is $\sum _{\left\{p,q\right\}\epsilon N}{u}_{\left\{p,q\right\}}·T\left(f\left(p\right)\ne f\left(q\right)\right),$

[0208] which is the cost of the multiway cut corresponding to f.

[0209] This shows that minimizing the Potts model energy E_P(f) on an arbitrary P and N is intractable. In computer vision, however, P is usually a planar grid, and combinatorial problems that are intractable on arbitrary graphs sometimes become tractable on the plane or grid.

[0210] The energy minimization problem is intractable even when restricted to a planar grid. The reduction is from a special case of the multiway cut problem, where G is a planar graph with degree 11 and all the edges have weight 1, which is NP-hard. G must be embedded in a grid of pixels, which happens in two stages. In the first stage, G is converted into a planar graph of degree 4. In the second stage, this graph is embedded in the grid. This embedding can be done in polynomial time; after it is done, each vertex v E G corresponds to a connected set of pixels S(v) in the grid, and the adjacency relationships among vertices in G has been preserved.

[0211] The proof now proceeds along the same lines as theorem 7.3, except for three subtleties. First, for every vertex v all pixels in S(v) are given the same label. This is ensured by making the edge weights K between adjacent pixels in S(v). Second, there are gaps when G is embedded in the grid. This can be solved by adding additional “grid pixels”, which D forces to have the extra label 0 (D will prevent non-grid pixels from having label 0 by making D_p(0)=K). The edge weights are taken between grid pixels and non-grid pixels to be one. The cost of a feasible labeling is now the cost of the corresponding multiway cut plus a constant. Third, the constant K>E_P(f*) must be now chosen more carefully.

[0212] Experimental Results

[0213] The algorithms of the present invention were applied to image restoration and visual correspondence and achieved promising results.

[0214] The experimental results are from visual correspondence for stereo and motion. In visual correspondence, two images are taken at the same time from different view points for stereo, and are also taken from the same view point but at different times for motion. For each pixel in the first image there is a corresponding pixel in the second image which is a projection along the line of sight of the same real world scene element. The difference in the coordinates of the corresponding points is called the disparity. In stereo the disparity is usually one-dimensional because corresponding points lie along epipolar lines. In motion the disparity is usually a two-dimensional quantity. The disparity varies smoothly everywhere except at object boundaries, and corresponding points are expected to have similar intensities. The correspondence problem can be formulated as the energy minimization problem $E\left(f\right)=\sum _{\left\{p,q\right\}\epsilon N}{V}_{p,q}\left(f_p,f_q\right)+\sum _{p\epsilon P}D\left(I_p-I_{p+f_p}^{\prime }\right)$

[0215] Here P is the set of all pixels in the first image, I_pis the intensity of pixel p in the first image, I_q′ is the intensity of pixel q in the second image, and p+f_pstands for the pixel with coordinates of p shifted by disparity f_p. The data penalty D is small if there is a small difference between I_pand i′_p+f.

[0216] For the experiments, three energy functions were used. The first energy function, called E_Q, uses the truncated quadratic V(f_p,f_q)=min(K,|f_p,f_q|²) (for some constant K) as its smoothness term. This choice of V does not obey the triangle inequality, so E_Qwas minimized using the swap algorithm of the present invention.

[0217] The second (E_P) and the third (E_L) energy functions use, correspondingly, the Potts model and the truncated L₂distance as their smoothness penalty V. Both of these obey the triangle inequality and E_Pand E_Lwere minimized with the expansion algorithm of the present invention.

[0218] Data Term

[0219] If the pixels p and q correspond, they are assumed to have similar intensities I_pand I′_q. Thus (I_p−I′_q)²is frequently used as a penalty for deciding that p and q correspond. This penalty has a heavy weight unless I_p≈I′_q. However there are special circumstances when corresponding pixels have very different intensities due to the effects of image sampling. Suppose that the true disparity is not an integer. If a pixel overlaps a scene patch with high intensity gradient, then the corresponding pixels may have significantly different intensities.

[0220] For stereo, a technique is used to develop a D_pthat is insensitive to image sampling. First it was measured how well p fits into the real valued range of disparities $\left(d-\frac{1}{2},d+\frac{1}{2}\right)$

[0221] by $C_{\mathrm{fwd}}\left(p,d\right)=\underset{d-\frac{1}{2}\le \times \le d+\frac{1}{2}}{\min }\left|I_{\alpha }-I_{p+d}^{\prime }\right|$

[0222] Fractional values I′_p+xare the result by linear interpolation between discrete pixel values. For symmetry, $C_{\mathrm{rev}}\left(p,d\right)=\underset{p-\frac{1}{2}\le \times \le p+\frac{1}{2}}{\min }\left|I_x-I_{p+d}^{\prime }\right|$

[0223] is also measured. C_fwd(p,d) and C_rev(p,d) can be computed with just a few comparisons. The final measure is

[0224] C(p,d)=(min{C_fwd(p,d),C_rev(p,d),Const})²

[0225] Here Const is used to make the measure more robust.

[0226] Static Cues

[0227] In this section, how to choose different V_p,qfor each pair of interacting pixels {p,q} to take advantage of contextual information is described. For simplicity, consider the case of the Potts model, i.e. V_p,q=u_{p,q}·T(f_p≠f_q). The intensities of pixels in the first image contain information that can significantly influence the assessment of disparities without even considering the second image. For example, two neighboring pixels p and q are much more likely to have the same disparity if it is known that I(p)≈I(q). Most methods for computing correspondence do not make use of this kind of contextual information.

[0228] Contextual information can be incorporated into the framework by allowing u_{p,q} to vary depending on their intensities I_pand I_q. Let

[0229] u_{p,q}=U(|I_p−I_q|)

[0230] Each u_{p,q} represents a penalty for assigning different disparities to neighboring pixels p and q. The value of the penalty u_{p,q} should be smaller for pairs {p,q} with larger intensity differences |I_p−I_q|. In practice, an empirically selected decreasing function U(•) is used. Note that instead of (19) the coefficients u_{p,q} could also be set according to an output of an edge detector on the first image. For example, u_{p,q} can be made small for pairs {p,q} where an intensity edge was detected and large otherwise. Segmentation results can also be used.

[0231] The following example shows the importance of contextual information. Consider the pair of synthetic images below, with a uniformly white rectangle in front of a black background.;

[0232] There is a one pixel horizontal shift in the location of the rectangle, and there is no noise. Without noise, the problem of estimating f is reduced to minimizing the smoothness term E_smooth(f) under the constraint that pixel p can be assigned disparity d only if I_p=I′_p+d.

[0233] If u_{p,q} is the same for all pairs of neighbors {p,q} then E_smooth(f) is minimized at one of the labeling shown in the picture below. Exactly which labeling minimizes E_smooth(f) depends on the relationship between the height of the square and the height of the background.

[0234] Suppose now that the penalty u_{p,q} is much smaller if I_p≠I_qthan it is if I_p=I_q. In this case the minimum of E_smooth(f) is achieved at the disparity configuration shown in the picture below. This result is much closer to human perception.

[0235] Real Stereo Imagery With Ground Truth

[0236] FIG. 10 shows the results from a real stereo pair with known ground truth. The left image of the real stereo pair is shown in FIG. 10(a). FIG. 10(b) shows the ground truth for this stereo pair.

[0237] For this stereo pair, E_Pwas used. The results were compared against annealing and normalized correlation. For normalized correlation, parameters which give the best statistics were chosen. Several different annealing variants were implemented. The one that gave the best performance was used. This was the Metropolis sampler with a linearly decreasing temperature schedule. To give it a good starting point, simulated annealing was initialized with the results from normalized correlation. In contrast, for the algorithms of the present invention, the starting point is unimportant. The results differ by less than 1% of image pixels from any starting point that was tried.

[0238] FIGS. 10(c), and (d) show the results of the swap and expansion algorithms for λ=20. FIGS. 10(e) and (f) show the results of normalized correlation and simulated annealing.

[0239] The table below summarizes the errors made by the algorithms. In approximately 20 minutes, simulated annealing reduces the total errors normalized correlation makes by about one fifth and it cuts the number of ±1 errors in half. It makes very little additional progress in the rest of 19 hours that it was run. The expansion, swap, and jump algorithms of the present invention make approximately 3 times fewer ±1 errors and approximately 5 times fewer total errors compared to normalized correlation.

[0240] The expansion and swap algorithms of the present invention perform similarly to each other. The observed slight difference in errors is quite insignificant (less than one percent). At each cycle the order of labels to iterate over is chosen randomly. Another run of the algorithms might give slightly different results, where expansion algorithm might do better than the swap algorithm. In general, very slight variation between different runs of an algorithm was observed. However the difference in the running time is significant. On average the expansion algorithm converges 3 times as rapidly as the swap algorithm.

[0241] FIG. 11 shows the graph of E_smoothversus time for the algorithms of the present invention versus simulated annealing. Notice that the time axis is on the logarithmic scale. The graph for E_datais not shown because the difference in the E_dataterm among all algorithms is insignificant, as expected from the following argument. Most pixels in real images have nearby pixels with very similar intensities. Thus for most pixels p there are a few disparities d for which D_p(d) is approximately the same and small. For the rest of d's, D_p(d) is quite large. This latter group of disparities is essentially excluded from consideration by energy minimizing algorithms. The remaining choices of d are more or less equally likely.

[0242] Thus the E_dataterm of the energy function has very similar values for the methods of the present invention and simulated annealing he methods of the present invention quickly reduce the smoothness energy to around 16,000, while the best simulated annealing can produce in 19 hours is around 30,000, nearly twice as bad.

[0243] The expansion algorithm gives a convergence curve significantly steeper than the other curves. In fact the expansion algorithm makes 99% of the progress in the first iteration.

[0244] The algorithms appear to be quite stable in the choice of parameter λ. For example the table in FIG. 12 gives the errors made by the expansion algorithm for different choices of λ. For small λ the algorithm makes a lot of errors because it overemphasizes the data, for large values of λ the algorithm makes a lot of errors because it overemphasizes the prior. However for a large interval of λ values the results are good.

[0245] SRI Tree Stereo Pair

[0246] In the well-known SRI stereo pair whose left image is shown in FIG. 13(a) the number of disparities is larger, and E_Pdoes not work as well. FIG. 13(b) and (c) compares the results of minimizing E_Pand E_L. Notice that there are fewer disparities found in FIG. 13(b), since the piecewise constant prior tends to produce large regions with the same disparity.

[0247] Motion

[0248] FIG. 14(a) shows one image of a motion sequence where a cat moves against moving background. This is a difficult sequence because the cat's motion is non-rigid. FIGS. 14(b) and (c) show the horizontal and vertical motions detected with the swap algorithm of the present invention. Notice that the cat has been accurately localized. Even the tail and parts of the legs are clearly separated from the background motion.

[0249] FIGS. 15aand 15bshows the output from the well-known flower garden sequence. Since the camera motion is nearly horizontal, the camera motion was simply delayed. E_Pwas used to achieve these results.

[0250] Bayesian Labeling of Markov Random Fields

[0251] In this appendix it is shown that minimization of the energy function in (1) is equivalent to computing the maximum a posteriori estimate of a Markov Random Field (MRF).

[0252] Let P be a set of sites, La set of labels, and N={{p,q}|p,q∈P} a neighborhood system on P. Let F=F₁, . . . , F_nbe a set of random variables defined on P. Each F_ptakes values in the label set L. A particular realization of the field will be denoted by f={f_p|p∈P}, which is also called a configuration of the field F. As usual, P(F_p=f_p) will be abbreviated by P(f_p). F is said to be a Markov Random Field (MRF) if:

[0253] (i) P(f)>0 ∀f∈F.

[0254] (ii) P(f_{p |fP−{p}})=P(f_p|f_N)

[0255] where P−{p} denotes set difference, f_Ndenotes all labels of sites in N_p, and F denotes the set of all possible labelings.

[0256] The easiest way to specify an MRF is by joint distribution. This theorem proves the equivalence between MRFs and Gibbs random fields.

[0257] Before defining Gibbs random fields, a clique needs to be defined. A set of sites is called a clique if each member of the set is a neighbor of all the other members. A Gibbs random field can be specified by the Gibbs distribution: $P\left(f\right)=Z^{-1}·\exp \left(-\sum _{c\epsilon C}{V}_c\left(f\right)\right),$

[0258] where C is the set of all cliques, Z is the normalizing constant, and {V_c(f)} are functions from a labeling to non-negative reals, called the clique potential functions.

[0259] Thus to specify an MRF, the clique potential functions need to be specified. A first order MRF is considered, which means that for all cliques of size larger than two the potential functions are zero, and for the cliques of size two the potential functions are specified by

[0260] V_c(f)=V_p,q(f_p,f_q).

[0261] This defines an MRF with the joint distribution: $P\left(f\right)=Z^{-1}·\exp \left(-\sum _{\left\{p,q\right\}\epsilon N}{V}_{p,q}\left(f_p,f_q\right)\right)$

[0262] In general, the field F is not directly observable in the experiment. A popular way to estimate its realized configuration {circumflex over (f)} based on an observation d is the maximum a posteriori (MAP) estimation. Using Bayes rule, the posterior probability can be written as $p\left(f❘d\right)=\frac{p\left(d❘f\right)p\left(f\right)}{p\left(d\right)}$

[0263] Thus the MAP estimate f* is equal to $\arg \underset{f\epsilon F}{\max }p\left(d❘f\right)p\left(f\right)=\arg \underset{f\epsilon F}{\min }\left(-\log p\left(d❘f\right)p\left(f\right)\right)$

[0264] Assume that he observation d_pat each pixel is independent and that

[0265] p(d_p|l)=C_p·exp(−D_p(l)) for l∈L,

[0266] where C_pis the normalizing constant, and D_pwas defined in Section 1. Then the likelihood can be written as $p\left(d❘f\right)\mathrm{\alpha exp}\left(-\sum _{p\epsilon P}{D}_p\left(f_p\right)\right)$

[0267] Writing out p(d) and p(d|f) with the above assumptions results in $f^{*}=\arg \underset{f\epsilon F}{\max }\sum _{\left\{p,q\right\}\epsilon N}{V}_{p,q}\left(f_p,f_q\right)+\sum _{p\epsilon P}{D}_p\left(f_p\right)$

[0268] which is the general form of the energy function being minimized.

[0269] It is to be understood that the above-described embodiments are simply illustrative of the principles of the invention. Various and other modifications and changes may be made by those skilled in the art which will embody the principles of the invention and fall within the spirit and scope thereof.