Superdifferential Cuts for Binary Energies
Tatsunori Taniai
University of Tokyo, Japan

Yasuyuki Matsushita
Osaka University, Japan

Takeshi Naemura
University of Tokyo, Japan

taniai@iis.u-tokyo.ac.jp

yasumat@ist.osaka-u.ac.jp

naemura@nae-lab.org

Abstract
We propose an efficient and general purpose energy
optimization method for binary variable energies usedin various low-level vision tasks. Our method can beused for broad classes of higher-order and pairwise nonsubmodular functions. We first revisit a submodularsupermodular procedure (SSP) [19], which is previously
studied for higher-order energy optimization. We thenpresent our method as generalization of SSP, which isfurther shown to generalize several state-of-the-art techniques for higher-order and pairwise non-submodular functions [2, 9, 25]. In the experiments, we apply our method toimage segmentation, deconvolution, and binarization, andshow improvements over state-of-the-art methods.

Initialization

pPBC [25]

FTR [10]Figure 1. Matching foreground color distribution using the proposed SDC-GEO, parametric pseudo bound cuts (pPBC) [25], andfast trust region (FTR) [10] with two types of initialization. pPBC
can only successively reduce the initial segment, while our method
allows arbitrary directions of optimization and is thus robust to initialization. (L2 distance for 643 bins of RGB histograms are used)ing energies by iteratively minimizing either first-order approximations (gradient descent approach) [10, 11] or upperbounds (bound optimization approach) [25, 2, 26, 20, 1]
of non-linear functions using graph cuts. The bound optimization approach has some advantages over the gradient descent approach [2]: It requires no parameters (e.g.,
step-size) and never worsens the solutions during iterations.
But we must in turn derive appropriate bounds for individual functions. A notable work is auxiliary cuts (AC) [2]by Ayed et al., where they derive general bounds for broadclasses of non-linear functionals for segmentation. However, the bounds derived in [2] are formulated to successively reduce target regions; thus the resulting segments are
restricted within initial segments. Such a property actuallylimits the applications and accuracy of the method.
In order to derive more accurate and useful bounds, werevisit a submodular-supermodular procedure (SSP) [19], ageneral bound optimization scheme for supermodular functions. We then propose a bound optimization method as
generalization of SSP. Unlike SSP, our method can be usedeven for non-supermodular functions; and unlike AC, it allows bi-directional optimization (see Fig. 1 for an illustration in segmentation) and can produce more accurate approximation bounds. We further show that our method canbe seen as generalization of AC and some state-of-the-art
method [9] for pairwise non-submodular functions.

1. Introduction
Many low-level vision problems such as image segmentation, binarization, denoising, and tracking are often formulated as binary energy minimization [4, 22, 1, 25, 6]. For
example, in image segmentation, the use of Markov random
field formulations [8] and graph cuts [17, 4] has been becoming one of primary approaches [3, 23, 21, 25, 9, 2, 10,
11, 26, 20, 1, 22]. In this approach, the energy function istypically formulated as
E(S) = R(S) + Q(S),SDC-GEO(1)

where R(S) describes appearance consistencies between
resulting segments S and given information about target regions, and Q(S) enforces smoothness on segmentboundaries. The form of R(S) is often restricted to simple linear (i.e., pixelwise unary) forms [3, 23, 21] becausegraph cuts allow globally optimal inference only for unaryand submodular pairwise forms of energies [17]. However,recent studies [25, 2, 10, 11, 26, 20, 1, 22] have shown thatthe use of higher-order information (or non-linear terms)
can yield outstanding improvements over conventional pixelwise consistency measures.
In general, higher-order terms involve difficult optimization problems. Recent promising approaches try reduc1

This paper makes the following contributions:• we propose an optimization method for broad classesof higher-order and pairwise non-submodular functions that allows arbitrary directions of convergenceand outperforms the state-of-the-art [9, 25, 10, 2].
• our method generalizes previous optimization methods
including from early [19] to state-of-the-art [9, 25, 2,
1] methods.

1.1. Scope of the Problems
Our objective is to seek the binary variable S such that it
minimizes E(S) of Eq. (1). In this paper, we focus on three
types of energy functions. Before defining those functions,
we define the following function as the basis of all types.
Pairwise Submodular Functions.
Let si ∈ {0, 1} be a binary variable defined for pixels i ∈ Ωin the image domain Ω, and we define S = {i|si = 1} asa segment in the domain Ω. If R(S) in Eq. (1) is a linear
product of a function hi = h(i) : Ω → R
X
R(S) = hh, Si =
hi si
(2)
i∈Ω

and Q(S) is the sum of pairwise functions
X
Q(S) =
mij si sj ,Type-2: Fractional Higer-Order Energies.
We further deal with non-supermodular functions:
X
Rtype2 (S) =fz (hgz , Si/hwz , Si)(6)zsuch that Rz (S)=fz (hgz , Si/hwz , S0 i) becomes a Type-1 term for fixed S0 .
Examplesof such Pfunctions include  the KL-divergence
R(S) = − z pz log nSz /|S| +  + const. and the
P p
Bhattacharyya coefficient R(S) = − z pz nSz /|S|,both
P are used for image segmentation [1, 2, 25]. Here,z pz = 1 is the target distribution.
Type-3: Pairwise Non-submodular Energies.
We also consider pairwise non-submodular energies, i.e.,
Q(S) of Eq. (3) containing non-submodular or supermodular terms (i.e. mij > 0). QPBO [16] is often used for suchfunctions, but it leaves many variables unlabeled when theamount of non-submodular terms is significant [9].

(3)2. Submodular-Supermodular Procedure(i,j)∈Ω

and if all the quadratic terms are non-positive (mij ≤ 0),then E(S) is submodular and can be globally minimized
via graph cuts [4, 17] in polynomial time.
Type-1: Higer-Order Supermodular Energies.
We consider the energies E(S) with pairwise submodular
Q(S) and the following form of R(S):X
X
Rtype1 (S) =Rz (S) =fz (hgz , Si)(4)z18] and tracking [13], where z ∈ Z is a bin of color orfeature histograms, hz is the given target histogram, and nSz
is the number of pixels in S that fall into the bin z. A volp
umetric constraint R(S) = |V0 − |S|| has been also used
for medical image analysis [11, 10].

z

Before presenting our method, we review SSP [19], an
optimization method for general supermodular functions,and later propose our method as its generalization.
SSP is classified as a bound optimization approach,where a tight upper bound function Ê(S|S t ) given an auxiliary variable S t is derived for E(S), i.e.,E(S) ≤ Ê(S|S t ) and E(S t ) = Ê(S t |S t ).

(7)Then, the bound is iteratively minimized as
S t+1 = arg min Ê(S|S t ), (t = 0, 1, 2, ...).

(8)+

where fz (x) is convex and gz (i) : Ω → R is a positive
function. Hence, Rz (S) is supermodular, i.e., if it satisfies
the following inequality for any X, Y ⊆ Ω:
Rz (X) + Rz (Y ) ≤ Rz (X ∩ Y ) + Rz (X ∪ Y ).

(5)

Note that Rz0 (S) = −Rz (S) is submodular, if Rz (S) is supermodular. While any submodular functions can be minimized in polynomial time [24], the minimization of supermodular functions is NP-hard.
Higher-order supermodular functions have been used,
for example,
as a Lp -distance histogram constraintPp
R(S) = z∈Z hz − nSz
for co-segmentation [22, 28,Here, it is guaranteed that the energy does not go
up, i.e. E(S t ) ≥ E(S t+1 ) holds for any t, because
E(S t ) = Ê(S t |S t ) ≥ Ê(S t+1 |S t ) ≥ E(S t+1 ).

2.1. Permutation Bounds
SSP derives tight bounds for supermodular functions
based on the superdifferential [7, 12], which is a similar
concept to the subderivative of continuous functions. Givena supermodular functions R(S) and S t ⊆ Ω, a modular (orlinear) function H(S) := hh, Si + R(∅) that satisfies
H(S) − H(S t ) ≥ R(S) − R(S t )

(9)

��

��
�(1)

�(1)

�1�

�(2)

�2�

�(3)

�(4)

�(5)

�(6)

�(1)

�(2)

�1�

�3�

�(2)

�(3)

�(3)

�(4)

�(4)
�(5)

�(6)

�2�

�3�
(b) Our generalized chain(a) Chain in SSP [19]Figure 2. Illustration of the chain and permutation.

�Segments
Geodesic [5]Euclidean
Figure 3. Examples of Euclidean and geodesic
distance [5] for segments (pink).

Unary cost ℎ� (�)ℎ(� 1 )

��ℎ(� 2 )

∅ �1� �2� �3� �4� �5� Ω��∅

�1� �2� �3�

+ �-term

Energy Value

Bound for �(�)Energy Value

Supermodular �(�)Energy Value

∞

Ω

∅

��

Ω

(b) Our generalized bound(a) Bound in SSP [19](c) Bound in AC [2]Figure 4. Illustration of upper bounds for supermodular functions. The supermodular function R(S) and its bounds are visualized by blue
and red lines, respectively. The green arrows in (a) show the unary costs h(σ(j)) of a supergradient H σ (S|S t ).

for any S ⊆ Ω, is called a supergradient of R at S t .
We denote ∂R(S t ) the set of all the supergradients of R
at S t , which is called the superdifferential. Notice that
if H(S t ) = R(S t ) holds, then H(S) gives a tight upperbound to R(S). Such extreme points of ∂R(S t ) may beobtained using the following theorem.
Theorem 1 (Theorem 6.11 in [7]):For any S t ⊆ Ω,a modular function H(S) = hh, Si + R(∅) is an extreme
point of ∂R(S t ), if and only if there exists a maximal chain
C : ∅ = S0 ⊂ S1 ⊂ · · · ⊂ Sn = Ω,

(10)

3. Proposed Methodwith Sj = S t for some j, such that
H(Sj )−H(Sj−1 ) = R(Sj )−R(Sj−1 ), (j = 1, ..., |Ω|).
(11)
Based on this theorem, SSP [19] derives a supergradientH σ (S|S t ) of R(S) at S t by the following greedy algorithm.
Let σ(j) : {1, 2, ..., |Ω|} → Ω be a permutation of Ω thatassigns the elements in S t to the first |S t | positions, i.e.,σ(j) ∈ S t if and only if j ≤ |S t |. A maximal chain C σ isthen defined as S0σ = ∅ and Sjσ = {σ(1), σ(2), ..., σ(j)},σtso S|St | = S . See Fig. 2 (a) for an illustration. Using thisσchain C , a supergradient H σ (S|S t ) is obtained as
H σ (S|S t ) = hhσ , Si + R(∅),

Figure 4 (a) illustrates how hσ (σ(j)) are computed. Since
R(S) is supermodular, variables in earlier positions of σ areassigned lower unary costs, i.e., more prone to be labeled
as 1 via cost-minimization. Therefore, if we knew a priori
how likely each variable is 1, the ideal permutation would
arrange variables in order of decreasing likelihood so as tomaximize the likelihood via cost minimization. Here, thebounds H σ (S) approximate R(S) tightly at along the chainof solutions {Sjσ }. However, because H σ (S) can largely
deviate from R(S) at other than {Sjσ }, SSP is problematicwhen likelihood or permutation is given inaccurate.

(12)

where each unary cost hσ (i) is given byσhσ (σ(j)) = R(Sjσ ) − R(Sj−1), (j = 1, 2, .., |Ω|). (13)In the following sections, we first present our key idea
by extending SSP [19]. We then show how to apply it toand optimize the three types of functions in Sec. 3.2, anddescribe implementation details in Sec. 3.3.

3.1. Grouped Permutation Bounds
In this section, we derive general bounds for Type-1
terms R(S) = f (hg, Si) by extending the SSP’s permutation scheme. First, we introduce a grouped permutation π, which is made by grouping SSP’s ordered-elements
σ = {σ(1), σ(2), · · · , σ(|Ω|)} into M (M ≤ |Ω|) groups:
π(1), π(2), · · · , π(M ) ⊆ Ω. Each group π(j) contains
some consecutive elements of σ: π(j) = {σ(j 0 ), σ(j 0 +1), · · · , σ(j 0 + m)}, and groups are mutually disjoint:
π(j) ∩ π(j 0 ) = ∅ if j 6= j 0 . Using this grouped permutation π we define a chain C π : S0π = ∅ and Sjπ =π(1) ∪ π(2) ∪ ... ∪ π(j) as illustrated in Fig. 2 (b). Here,we make sure that any group does not across σ(|S t |) andσ(|S t | + 1), so there exists Sjπ = S t for some j. Then,our bound for R(S) is defined similarly to that of SSP in
Eq. (12) as
H π (S|S t ) = hhπ , Si + R(∅),(14)

where unary costs hπ (i) : Ω → R are defined for i ∈ π(j)and j = 1, 2, ..., M asπhπ (i) = g(i) R(Sjπ ) − R(Sj−1) / hg, π(j)i .
(15)
The essence of this formulation becomes clearer if we assume g = 1 so that H π (S|S t ) becomes piecewise-meanapproximations of SSP’s bound H σ (S|S t ), which is visualized in Fig. 4 (b) using the example permutation π
shown in Fig. 2 (b). Note that our bound H π (S|S t ) becomes equivalent to the supergradients H σ (S|S t ) of SSP, if
π(j) = {σ(j)}.
Proposition 1: The function H π (S|S t ) satisfies the conditions of Eq. (7) and is thus a tight upper bound for R(S).
Proof. See our supplementary.
The spirit behind this grouping or piecewise-meanapproximation scheme is to make fine/coarse approximation bounds when permutations σ are accurate/inaccurate.
Specifically, when the permutation of σ(j) and σ(j + 1)
is unreliable, we put them into the same group in order totreat them equally and leave a decision (i.e., which is more
likely to be labeled as 1) to other interactions, e.g., pairwisesmoothness terms. As we will show in Sec. 3.3, our method
makes coarse-to-fine approximation bounds by iterations.

3.2. Optimization Procedure
We optimize E(S) by iteratively minimizing its approximation function Ê(S|S t−1 ) derived by our grouped permutation bounds. Here, the minimization of Ê(S) is achieved
using a max-flow/min-cut algorithm [4].
Type-1: We derive a bound HzπP
(S|S t ) for each Rz (S)of Rtype1 (S), and set Ê(S|S t ) = z Hzπ (S|S t ) + Q(S).
Here, it is guaranteed that minimization of Ê(S) does notincrease the energy E(S). Therefore, its optimization procedure is a simple iteration algorithm shown in Algorithm 1(without lines 5 and 6).
Type-2: Similarly to [2, 25], we approximate Rtype2 (S)by partially fixing S at S t as
XR̃type2 (S|S t ) =fz hgz , Si/ wz , S t .
(16)zHere, R̃type2 (S|S t ) is a Type-1 function. Therefore,
we can approximate R̃type2 (S|S t ) using our bounds

Algorithm 1 O PTIMIZATION F OR T YPE -1,3 [T YPE -2]
1: Initialize S 0
2: for t = 0, 1, 2, · · · do
3:
Create a permutation π // for all Types
4:
S t+1 ← argmin Ê(S|S t ) // for Type-1,35:[{S λ } ← argmin Ê(S|S t ) + λ(|S t − |S|) for ∀ λ]6:
[S t+1 ← argmin E(S) for S ∈ {S λ , S t }] // for Type-2
7: end for

H π (S|S t ). As long as S ⊆ S t is forced, this linear function H π makes a tight bound for Rtype2 [2].
However, we rather do not restrict S for allowing bidirectional optimization. In this case, the minimization ofÊ(S|S t ) = H π (S|S t ) + Q(S) may increase the actual energy E(S). For this we use the pseudo bound optimization
scheme of [25], where we make a family of relaxed boundsÊλ (S) = Ê(S|S t ) + λ(|S t | − |S|), and we exhaustively
search S λ = argmin Êλ (S) for all λ ∈ (−∞, ∞) by usingparametric maxflow [15]. Then we choose the best S λ thatminimizes E(S). Therefore, the optimization procedure in
Algorithm 1 takes lines 5 and 6 instead of 4 for Type-2.
Type-3: We make a bound Ê(S) by approximating eachof non-submodular terms R(S) = mij si sj (mij > 0) witha linear function hi si + hj sj + const. Its form depends on
both current values (sti , stj ) and the permutation π of {i, j},
and given as h(i0 ) = R(Sjπ0 ) − R(Sjπ0 −1 ) /|π(j 0 )| similarly to Eq. (15). We summarize the conversion in Table 1.

3.3. Implementation Details
We make permutations σ and π according to the signed
distance from the boundary of S t . In [22], the Euclidean
distance was used, but we use geodesic distance [5] formore effectively creating the permutations. Figure 3 shows
examples of both distances. We denote the geodesic distance for segments S by D(i|S, I) : Ω → R, and D(i|S) ≤0 for i ∈ S. See Eq. (7) of [5] and also our supplementary
for its definition. As discussed in [5], D(i|S, I) is efficiently
computed in O(|Ω|) using an approximate algorithm [27].
Using the geodesic distance, we construct the boundH π (S|S t ) in each iteration as follows. Firstly, we compute
D(i|S t ) for the current segments S t . Secondly, we make apermutation σ such that D(σ(j)) ≤ D(σ(j + 1)). Finally,
we make a grouped permutation π from σ. We process σ(j)from σ(2) to σ(|Ω|), and put σ(j) into the same group withσ(j −1) if D(σ(j))−D(σ(j − 1)) ≤ τ , while making surethe group does not across σ(|S t |) and σ(|S t |+1). Basically,
the size of the threshold τ reflects how much the permutation σ by D(i) is reliable. We empirically use a grouping
threshold given by τ = (µ + as)/(t + 1)κ (t = 0, 1, 2, ...),
where µ and s are the mean and standard deviation of distance differences |D(σ(j)) − D(σ(j − 1))|. We use thismonotonically decreasing thresholds, because as iterations

Table 1. The proposed linear conversions of non-submodular terms mij si sj with comparisons to the conversions of SSP [19] and LSA [9].
Below, the subscript of mij is omitted. Notice that our linear conversion can bee seen as generalization of both SSP and LSA-AUX.

(sti , stj )Ours(0, 0)mm
2 si + 2 sj if π(1) = {i, j}
msj if π(1) = {i}(0, 1)(1, 0)(1, 1)

msi if π(1) = {j}
msi
msj
m
m
s
+
s
if π(1) = {i, j}i
j22
msj if π(1) = {i}msi if π(1) = {j}SSP [19]

LSA-AUX [9]LSA-TR [9]msj if σ(1) = i
msi if σ(1) = jmm
2 si + 2 sj

0

msi
msj

msi
msj

msi
msj

msj if σ(1) = i
msi if σ(1) = jmm
2 si + 2 sj

msi + msj − mproceed the segments S t are expected to be more accurate
and so permutations σ by D(i|S t ) becomes accordingly
more reasonable. Also, when S 0 = Ω or S 0 = ∅, so
D(i|S t ) cannot be defined, we set π(1) = Ω. This makes
the full linear approximations of R(S) drawn from R(∅) toR(Ω), which give reasonable initial approximations.

4. Relationship with Prior Art
In this section, we discuss relationships between our
method and other previous methods [2, 25, 9].

4.1. Auxiliary Cuts (AC) [2] and pPBC [25]The authors of [2] derived a bound for Rtype1 (S) by using the Jensen’s inequality and assuming S ⊆ S t . Its essential form is given ash1, Si, (17)Aα (S|S t ) = Bshr (S|S t ) + αR(S t ) 1−
h1, S t i
where α ≥ 0 is a parameter and Bshr (S|S t ) is expressed as
Bshr (S|S t ) = R(∅) +R(S t ) − R(∅)
hg, Si .
hg, S t i

Hf ull (S|S t ) = Bshr (S ∩S t |S t )+Bexp (S ∩ S̄ t |S t ). (20)
Notice that it no longer requires the restrictions for S (i.e.,
S ⊆ S t or S t ⊆ S) used in [2, 25], which are turned out tobe unnecessary by our derivation. Also notice that Hf ulldoes not depend on the permutation π. Although Hf ulldoes not accurately approximate its original function R(S),it is still useful when permutations are inaccurate (e.g. atfirst iterations). Our method is designed to behave between
AC [2] and SSP [19] and to produce coarse-to-fine approximation bounds as iterations proceed, by using the monotonically decreasing grouping-threshold τ defined in Sec. 3.3.

4.2. Local Submodular Approximations (LSA) [9](18)

From Eq. (18), Bshr (S|S t ) can be seen as the linear approximations of R(S) drawn from R(∅) to R(S t ). See
Fig. 4 (c), where the solid red line visualizes this linear approximation bound, the infinite bound reflects the restriction
S ⊆ S t , and the dotted line depicts the effect of the α-termin Eq. (17). The use of this bound results in successively
shrinking segments: S 0 ⊇ S 1 ⊇ S 2 . pPBC [25] extends
AC by exhaustively searching the best α ∈ (−∞, ∞) ineach iteration using parametric maxflow [15].
To point out a relationship to our method, we use another bound Bexp (S|S t ) for R(S), which is similar to
Bshr (S|S t ) and can be derived by [2]’s derivations as
Bexp (S|S t ) = R(S t )+This bound can be seen as the linear approximation of R(S)drawn from R(S t ) to R(Ω).
We now show that these bounds Bshr (S|S t ) and
Bexp (S|S t ) can be derived using our grouped permutation
bounds. When a grouped permutation is given as π(1) = S tand π(2) = S̄ t with S̄ t := Ω \ S t , our bound H π (S|S t ) becomes the following form:

R(Ω) − R(S t )g, S \ S t . (19)hg, Ω \ S t i

Here, S t is restricted to S t ⊆ S, so the iterative minimization of Bshr (S|S t ) successively expands the segments S.

Very recently, an optimization method for pairwise nonsubmodular energies called LSA [9] has been proposed andshown to outperform other state-of-the-art methods suchas TRW-S [14] and QPBO [16]. This method approximates non-submodular pairwise terms by linear terms. In[9], two types of approximation conversions are proposed.
LSA-TR applies a Taylor-based approximation and uses thegradient-descent framework of FTR [10]. LSA-AUX uses abound-based approximation and uses the bound optimization framework. The linear conversions of LSA-TR andLSA-AUX are summarized in Table 1.
As you can see from the table, our linear conversion
includes the conversion of LSA-AUX. In fact, our conversion with a full-grouping π produces the same boundswith LSA-AUX. Furthermore, as mentioned in [9], there
are other types of bounds that can be made by the permutation scheme of SSP [19]. Our method generalizes bothconversions of SSP and LSA-AUX, and adaptively chooses
either conversion for each term.

5. Experiments
For evaluation, we use four variations of our method:
SDC-GEO is the proposed method described in Sec. 3.
SDC-DIST uses the standard Euclidean distance for making permutations σ, but the other settings are the samewith SDC-GEO. SSP-DIST is SSP [19] that follows theimplementations by Rother et al. [22]. It is basically thesame with SDC-DIST but uses no mean approximations forbound constructions. We also use SSP-GEO, which is thesame with SSP-DIST but uses the geodesic distance. Notethat when S 0 = Ω at the first iterations, we make permutations σ randomly for SSP-DIST and SSP-GEO, based on
10 × 10-pixels of patches as described in [22].
We also compare with three state-of-the-art methods
for higher-order energies: AC [2] uses the bound Aα of
Eq. (17). pPBC [25] uses the same bound but exhaustively chooses the best α using parametric maxflow [15].
FTR [10] is the state-of-the-art of gradient-descent methods. For pairwise energies, we compare with LSA-AUX,TR [9] and pPBC-T,-B,-L [25] as state-of-the-art.
All the methods are implemented by C++, and run on asystem with a 3.5GHz Core i7 CPU and 16GB RAM.

5.1. Segmentation via Distribution Matching
Similarly to [25, 2, 22], we evaluate the performancesof the methods using the GrabCut dataset [23]. For pureevaluations of optimization performances, we learn the target histograms from the ground truth1 . We use a standard 16-neighbor
pairwise smoothness term similar to [23]:P
Q(S) = λ ij max(wij , )|si − sj |/|pi − pj | where pi is
pixel coordinates and wij = exp(−β|Ii − Ij |2 ). Here, β is
automatically estimated as β = 2E[|Ii − Ij |2 ] the expectation over all neighbor pairs. We use RGB-histograms.
Type-1: L2 and L1 Distance of Histograms
For Type-1 terms, we use the L2 and L1 distances between
histograms. For {λ, }, we use {1.0, 0.5} and {0.5, 0.5}
for the L2 and L1 -distances, respectively. S t is trivially
initialized as S 0 ← Ω. For both SDC-GEO and SDC-DIST,
we use a grouping threshold τ = (µ + 8s)/(t + 1)2.5 .
Tables 2 and 3 summarize the performance comparisons,showing average misclassified pixel rates, energy values ofE(S) and R(S), running times, and individual-image comparisons with SDC-GEO. Among the seven methods, ourproposed method SDC-GEO outperforms the others for all
error and energy scores. In some cases, SDC-GEO completely outperforms pPBC, AC, FTR, and SSP-DIST for all
individual images. Note that the error rates of SSP-DIST
are better than the rates originally reported in [22] because
1 If the target histograms are inaccurate, the minimum solutions ofE(S) are deviated from the ground truth [26], and the error rate criteria
thus does not reflect the actual performances of the optimization methods.

Table 4. Evaluations on the GrabCut dataset [23] using the Bhattacharyya distance and KL divergence. We use 643 bins and thebounding box initialization.
Error (%)E(S)
Time (sec)MethodBhat. KL-div. Bhat. KL-div. Bhat. KL-div.
SDC
0.373 0.515 -14906 7025 19.8 12.9
pPBC [25]0.498 0.818 -14870 7057
1.8
1.5
AC [2]
18.29 16.07 -11878 7490
0.4
0.4
FTR [10]0.435 1.076 -14894 7049
2.9
6.4

the definition of pixel feature histograms is different2 . Comparing the results of SDC-DIST and SSP-DIST with L2 and643 bins, SDC-DIST finds more accurate segmentations inspite of the higher energies. This is because in SSP-DIST
the appearance consistencies are forced regardless of how
permutations σ and corresponding bounds are inaccurate,
resulting in highly non-smooth, visibly bad local minimas.
In Fig. 5, we show example results of L2 and L1 with 64
bins. Figures 6 and 7 show the plots of the accuracy transitions using the L2 and L1 -distances w.r.t. the number ofbins. As shown, SDC-GEO is robust to the difference of thenumber of bins. In order to show robustness to initialization,
we compare SDC-GEO, pPBC, and FTR using two types of
initialization shown in Fig. 1. Unlike pPBC (and AC) thatcan only reduce the target regions, our method is robust to
initialization and finds very accurate solutions even for suchdifficult camouflage images.
Type-2: Bhattacharyya Distance and KL Divergence
We use Bhattacharyya distance and KL divergence as
Type-2. As described in Sec. 3.2, we use the pseudo boundoptimization scheme of [25] using parametric maxflow [15]for our method SDC-GEO. We show the performance comparisons with pPBC, AC, and FTR in Table 4, where ourmethod outperforms the others in error and energy values. It
is worth noting that although SSP [19] is originally oriented
for supermodular functions, our extended method successfully optimizes non-supermodular functions (Type-2).

5.2. Type-3: Image Deconvolution
˜
We make a blurred image
P I by a mean filter and additive
Gaussian noises: I˜i = 19 j∈Wi Ij +N (0, σ 2 ), where Wi isa 3 × 3 window centered at i.PWe recover the
P original image
I by minimizing E(S) = i∈Ω (I˜i − 19 j∈Wi sj )2 . In
Fig. 8, we show example results of our method, pPBC [25],and LSA [9] for two noise levels σ = 0.15, 0.30.
P In Fig. 9,
we show the plots of energies, squared errors ( i |si −Ii |2 ),and running times. LSA-TR and pPBC-T,-L reach the lower
energies but inaccurate results. Our method and LSA-AUX
perform best in terms of both accuracy and efficiency.
2 [22] uses a normalized 2D color vector and texton as a pixel feature,since the method is intended for co-segmentation and image retrieval.

Table 2. Evaluations on the GrabCut dataset [23] using L2 -distance. We show average error rates, E(S), R(S), times over 50 images. The
last column shows the number of images for which the proposed method (SDC-GEO) outperforms each method. We use 1923 and 643
bins.
MethodError (%)E(S)
R(S)
Time (sec)
SDC-GEO vs(L2 -distance)
1923
643
1923
643
1923
643
1923
643
1923
643
ref. Ground Truth
0
0
3569
3569
0
0
SDC-GEO
0.095
0.288
3511
3809
179
235
6.2
11.7
SDC-DIST
0.134
1.123
4121
18922
577
9946
2.3
2.7
32
46
SSP-GEO
0.116
1.132
3594
6477
166
377
4.9
15.1
24
45
SSP-DIST [22]
0.312
2.575
4415
13969
295
877
2.0
4.2
32
50
pPBC [25]
1.062
2.677
19381
19033213187
176315
36.7
24.6
49
50
AC [2]
1.214
3.542
19888
19572913646
177947
1.0
1.0
50
50
FTR [10]
1.859
3.167
21003
153212
17669
146394
36.8
121.2
50
48
Table 3. Evaluations on the GrabCut dataset [23] using L1 -distance. See Table 2 for the descriptions.
MethodError (%)E(S)
R(S)
Time (sec)
SDC-GEO vs(L1 -distance)
1923
643
1923
643
1923
643
1923
643
1923
643
ref. Ground Truth
0
0
1785
1785
0
0
SDC-GEO
0.033
0.205
1804
1882
54
120
5.4
9.4
SDC-DIST
0.041
0.242
1813
1967
70
229
1.9
3.8
38
42
SSP-GEO
0.043
0.943
1813
2766
31
298
3.7
9.7
28
46
SSP-DIST [22]
0.075
1.341
1898
3666
50
411
1.6
2.4
34
48
pPBC [25]
0.154
0.583
2013
2632
309
997
17.4
18.9
48
48
AC [2]
0.339
1.281
2502
4336
762
2476
0.6
0.6
50
50
FTR [10]
0.147
0.366
1908
2105
277
495
46.1
100.7
46
41

Ground Truth

SDC-GEO

SSP-DIST [22]pPBC [25]AC [2]FTR [10]3

Figure 5. L2 (top) and L1 (bottom) examples. The segments are initialized as all foreground, and 64 bins of histograms are used.
5.0F

2.5F
SDC-GEO
pPBC
AC
FTR

3.0F

SDC-GEO
pPBC
AC
FTR

2.0F
ErrorCRateC(F)ErrorCRateC(F)

4.0F

2.0F

1.0F

1.5F
1.0F

0.5F

0.0F

0.0F

32

64
96
128
160
NumberCofCBinsCperCChannel

192

Figure 6. Error rate transitions w.r.t. the number of bins (L2 )

32

64
96
128
160
NumberCofCBinsCperCChannel

192

Figure 7. Error rate transitions w.r.t. the number of bins (L1 )InputpPBC-T [25]SDC-GEO

pPBC-B [25]pPBC-L [25]LSA-AUX [10]LSA-TR [10]Figure 8. Image deconvolution results for two images with noise levels of σ = 0.15 (top) and σ = 0.30 (bottom).

EnergyT(10e+2)

10
8
6
4

SDC-GEO

PBC-B

SquaredPErrorPR10e+2)

SDC-GEO
pPBC-T
pPBC-B
pPBC-L
LSA-AUXLSA-TR

LSA-AUX

2
0
0.05

0.1

0.15
0.2
NoiseTLevelTσ

0.25

0.3

10
9
8
7
6
5
4
3
2
1
0

4.5

SDC-GEO
pPBC-T
pPBC-B
pPBC-L
LSA-AUXLSA-TR

SDC-GEO
pPBC-T
pPBC-B
pPBC-L
LSA-AUXLSA-TR

4.0
3.5

Time (sec)12

SDC-GEO

PBC-B
LSA-AUX

3.0
2.5

2.0
1.5
1.0
0.5
0.0

0.05

0.1

0.15
0.2
NoisePLevelPσ

0.25

0.3

0.05

0.1

0.15
0.2
Noise Level σ

0.25

0.3

Input

SDC-GEO

EnergyWS10eX4A

Figure 9. Performance comparisons of image deconvolution. Convergence energies, square errors, and running times w.r.t. noise levels areshown. The values are averaged over 30 random noise images at each point. We use τ = (µ + s)/(t + 1)2.5 for our method.

LSA-AUX [9]pPBC-T [25]pPBC-B [25]pPBC-L [25]Figure 10. Results of curvature optimization at the weight of 25.

6.4
6.3
6.2
6.1
6.0
5.9
5.8
5.7
5.6
5.5
5.4

SDC-GEO
pPBC-T
pPBC-B
pPBC-L
LSA-AUXLSA-TR

19 20 21 22 23 24 25 26 27 28 29 30
CurvatureWRegularizationWWeight

Figure 11. Convergence energies w.r.t. curvature regularization weights. We use τ = µ/(t + 1)2.5 .

5.3. Type-3: Curvature Regularization
We apply our method to a curvature regularization modelof [6] in image binarization. We show the input image andresults by our method, pPBC [25], and LSA [9] in Fig. 10.
The plots in Fig. 11 show energies at convergence w.r.t. regularization weights. Our method always reaches the lowest
energies and is most stable among all methods.

6. Conclusions
In this paper we have revisited SSP [19], an early
approach to higher-order energy optimization, and pro-

posed our method as generalization of SSP. The key ideaof our method is piecewise mean approximation bounds,
which are designed to produce coarse-to-fine approximation
bounds during iterations. We further show that our method
has close connections to some state-of-the-art methods [2,25, 9]. Although the proposed method shows promising improvements over state-of-the-art methods [25, 10, 9, 2], wewould like to further push the envelope by improving thedefinition of geodesic distance and the thresholding scheme
for making grouped permutations.

Acknowledgments
The authors would thank to anonymous reviewers forproviding valuable feedbacks for improving this paper. This
work was supported by JSPS KAKENHI Grant Number
14J09001.

References[1] I. B. Ayed, H. M. Chen, K. Punithakumar, I. Ross, and S. Li.
Graph Cut Segmentation with a Global Constraint: Recovering Region Distribution via a Bound of the Bhattacharyya
Measure. In Proc. of IEEE Conf. on Computer Vision and
Pattern Recognition (CVPR), pages 3288–3295, 2010. 1, 2
[2] I. B. Ayed, L. Gorelick, and Y. Boykov. Auxiliary Cuts
for General Classes of Higher Order Functionals. In Proc.
of IEEE Conf. on Computer Vision and Pattern Recognition(CVPR), pages 1304–1311, 2013. 1, 2, 3, 4, 5, 6, 7, 8
[3] Y. Boykov and M. P. Jolly. Interactive Graph Cuts for Optimal Boundary & Region Segmentation of Objects in N-D
Images. In Proc. of Int’l Conf. on Computer Vision (ICCV),
pages 105–112, 2001. 1[4] Y. Boykov and V. Kolmogorov. An Experimental Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision. IEEE Trans. Pattern Anal. Mach. Intell.
(TPAMI), 26:1124–1137, 2004. 1, 2, 4
[5] A. Criminisi, T. Sharp, and A. Blake. GeoS: Geodesic Image Segmentation. In Proc. of European Conf. on Computer
Vision (ECCV), pages 99–112, 2008. 3, 4
[6] N. El-Zehiry and L. Grady. Fast Global Optimization of Curvature. In Proc. of IEEE Conf. on Computer Vision and Pattern Recognition (CVPR), pages 3257–3264, 2010. 1, 8
[7] S. Fujishige. Submodular Functions and Optimization, volume 58. Elsevier Science, 2005. 2, 3
[8] S. Geman and D. Geman. Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images. IEEE
Trans. Pattern Anal. Mach. Intell. (TPAMI), 6(6):721–741,
1984. 1[9] L. Gorelick, Y. Boykov, O. Veksler, I. Ben Ayed, and A. Delong. Submodularization for Binary Pairwise Energies. In
Proc. of IEEE Conf. on Computer Vision and Pattern Recognition (CVPR), 2014. 1, 2, 5, 6, 8
[10] L. Gorelick, F. R. Schmidt, and Y. Boykov. Fast Trust Region for Segmentation. In Proc. of IEEE Conf. on Computer
Vision and Pattern Recognition (CVPR), pages 1714–1721,
2013. 1, 2, 5, 6, 7, 8
[11] L. Gorelick, F. R. Schmidt, Y. Boykov, A. Delong, and A. D.
Ward. Segmentation with Non-Linear Regional Constraints
via Line-Search Cuts. In Proc. of European Conf. on Computer Vision (ECCV), pages 583–0597, 2012. 1, 2
[12] R. Iyer, S. Jegelka, and J. Bilmes. Fast Semidifferentialbased Submodular Function Optimization. Proc. of International Conference on Machine Learning (ICML), 28(3):855–
863, 2013. 2[13] H. Jiang. Linear Solution to Scale Invariant Global Figure
Ground Separation. In Proc. of IEEE Conf. on Computer Vision and Pattern Recognition (CVPR), pages 678–685, 2012.
2[14] V. Kolmogorov. Convergent Tree-Reweighted Message
Passing for Energy Minimization. IEEE Trans. Pattern Anal.
Mach. Intell. (TPAMI), 28(10):1568–1583, 2006. 5[15] V. Kolmogorov, Y. Boykov, and C. Rother. Applications of
Parametric Maxflow in Computer Vision. In Proc. of Int’l
Conf. on Computer Vision (ICCV), 2007. 4, 5, 6
[16] V. Kolmogorov and C. Rother. Minimizing Nonsubmodular
Functions with Graph Cuts – A Review. IEEE Trans. PatternAnal. Mach. Intell. (TPAMI), 29(7):1274–1279, 2007. 2, 5
[17] V. Kolmogorov and R. Zabin. What Energy Functions Can
Be Minimized via Graph Cuts? IEEE Trans. Pattern Anal.
Mach. Intell. (TPAMI), 26(2):147–159, 2004. 1, 2
[18] L. Mukherjee, V. Singh, and C. Dyer. Half-Integrality based
Algorithms for Cosegmentation of Images. In Proc. of IEEE
Conf. on Computer Vision and Pattern Recognition (CVPR),
pages 2028–2035, 2009. 2[19] M. Narasimhan and J. A. Bilmes.
A SubmodularSupermodular Procedure with Applications to Discriminative Structure Learning. In Proc. of Conf. on Uncertainty in
Artificial Intelligence (UAI), pages 404–412, 2005. 1, 2, 3,
5, 6, 8[20] V.-Q. Pham, K. Takahashi, and T. Naemura. ForegroundBackground Segmentation using Iterated Distribution
Matching. In Proc. of IEEE Conf. on Computer Vision and
Pattern Recognition (CVPR), pages 2113–2120, 2011. 1[21] B. L. Price, B. Morse, and S. Cohen. Geodesic Graph Cut
for Interactive Image Segmentation. In Proc. of IEEE Conf.
on Computer Vision and Pattern Recognition (CVPR), pages
3288–3295, 2010. 1[22] C. Rother, V. Kolmogorov, T. Minka, and A. Blake.
Cosegmentation of Image Pairs by Histogram Matching–
Incorporating a Global Constraint into MRFs. In Proc. ofIEEE Conf. on Computer Vision and Pattern Recognition(CVPR), pages 993–100, 2006. 1, 2, 4, 6, 7
[23] C. Rother, V. Kolmogorv, and A. Blake. GrabCut: Interactive
Foreground Extraction Using Iterated Graph Cuts. Proc. of
SIGGRAPH (ACM Trans. on Graph.), 23:309–314, 2004. 1,6, 7[24] A. Schrijver. A Combinatorial Algorithm Minimizing Submodular Functions in Strongly Polynomial Time. J. of Combinational Theory Ser. B, 80(2):346–355, 2000. 2[25] M. Tang, I. Ben Ayed, and Y. Boykov. Pseudo-Bound Optimization for Binary Energies. In Proc. of European Conf. on
Computer Vision (ECCV), pages 691–707, 2014. 1, 2, 4, 5,
6, 7, 8[26] T. Taniai, V.-Q. Pham, K. Takahashi, and T. Naemura. Image
Segmentation using Dual Distribution Matching. In Proc.
of British Machine Vision Conf. (BMVC), pages 74.1–74.11,
2012. 1, 6
[27] P. J. Toivanen. New Geodesic Distance Transforms for Grayscale Images. Pattern Recogn. Lett., 17(5):437–450, 1996. 4[28] S. Vicente, V. Kolmogorov, and C. Rother. Cosegmentation
Revisited: Models and Optimization. In Proc. of European
Conf. on Computer Vision (ECCV), pages 465–479, 2010. 2