# Discontinuity Preserving Image Registration through Motion Segmentation: A Primal-Dual Approach

Authors: Silja Kiriyanthan; Ketut Fundana; Tahir Majeed; Philippe C. Cattin
Journal: Computational and Mathematical Methods in Medicine
 (2016)

---

## Abstract

Image registration is a powerful tool in medical image analysis and facilitates the clinical routine in several aspects. There are many well established elastic registration methods, but none of them can so far preserve discontinuities in the displacement field. These discontinuities appear in particular at organ boundaries during the breathing induced organ motion. In this paper, we exploit the fact that motion segmentation could play a guiding role during discontinuity preserving registration. The motion segmentation is embedded in a continuous cut framework guaranteeing convexity for motion segmentation. Furthermore we show that a primal-dual method can be used to estimate a solution to this challenging variational problem. Experimental results are presented for MR images with apparent breathing induced sliding motion of the liver along the abdominal wall.


---

## Body

## 1. Introduction

Image registration became an indispensable tool for many medical applications including the Image-Guided Therapy systems. Today’s image registration methods [1] are powerful at handling rigid as well as nonrigid motion on mono- and multimodal images. With the introduction of imaging technologies capable of capturing 4D organ motion, such as 4D-MRI [2], 4D-CT, or 4D-Ultrasound, imagery of sliding organs became a significant focus of research. Most of the state-of-the-art image registration methods cannot yet properly deal with these data sets, as their regularization causes continuous motion fields over the organ boundaries.In recent publications, the concept of the direction-dependent regularization of the displacement field has been introduced; see, for example, [3]. The drawback of these methods is the necessity of providing a good manual segmentation of the boundaries. Approaches that do not rely on prior manual segmentations have been a topic of research for some decades in classical computer vision such as optical flow. Transfer of these research results into the medical field was however scarce. Already in 1989, Mumford and Shah [4] proposed in their pioneering work a functional for image segmentation that avoids spatial smoothing in certain locations of the image, thus preserving discontinuities. Some years later, Weickert and Schnörr [5] proposed an extension of nonquadratic variational regularization for discontinuities preserving optical flow that uses spatiotemporal regularizers instead of flow-driven spatial smoothness. The resulting convex function guaranteed global convergence that could be solved with standard gradient based optimisation schemes. Vese and Chan [6] then introduced a level set framework based approach to efficiently solve the Mumford and Shah minimization problem for segmentation. Another influential approach based on the Total Variation (TV) norm, known to preserve discontinuities, was proposed by Rudin et al. [7], the ROF model. The beneficial behaviour of the TV norm was also exploited in recent registration and optical flow methods, as, for example, by [8, 9].Motion segmentation has been a topic of research for quite some time. In 1993, Nesi [10] proposed a variational optical flow approach that incorporates a variable for the presence of discontinuities and leads to a piecewise smooth motion estimation. Cremers and Schnörr introduced in [11] a variational approach seamlessly integrating motion segmentation and shape regularization into one single energy functional. In [12], the authors formulated the motion segmentation problem in a Bayesian inference framework, whereas the motion discontinuity could be represented either as an explicit spline function or by an implicit level set formulation allowing for an arbitrary number of objects to segment. As motion discontinuities mainly appear at object boundaries, it seems natural to combine motion segmentation and registration. Given a good motion segmentation, a proper discontinuous motion field can be estimated and vice versa. Amiaz and Kiryati [13] tried to leverage this property by combining motion segmentation with the optical flow method of Brox et al. [8] in a level set framework. They showed that motion segmentation leads to better discontinuities in the motion field. However, their approach highly depends on the initialisation for the motion segmentation and displacement field. Furthermore, the level set approach is prone to local minima and a rather good initialisation has to be provided in advance.In this paper we propose an elastic registration approach that handles discontinuities in the motion field by combining motion segmentation and registration in a variational scheme. For motion segmentation we make use of the so-called “continuous cuts” [14, 15] scheme that guarantees a globally optimal solution for a fixed motion field. In contrast to our previous work [16], the proposed energy functional is this time TV-L
1 regularized for both the motion segmentation and the displacement field and the fidelity term is formulated as a sum of absolute values of the constraints. Minimizing the TV-L
1-regularized functionals is, however, inherently difficult due to the nonsmoothness of the TV. This is an active field of research, such as [17–20] and the references therein which are mainly based on operator splitting methods in convex analysis that split energy functionals into the nonlinear and linear terms. We show how to solve the proposed complex energy functional with the primal-dual method of Chambolle and Pock [17], which is very efficient for a wide class of nonsmooth problems. The proposed method is then applied to register 2D MR image pairs with apparent breathing induced sliding motion of the liver along the abdominal wall. The preliminary version of this paper has been published in [21].


## 2. Method

In this section we describe the proposed registration method which integrates the displacement field estimation into the convex segmentation method of Chan et al. [14], in order to find smooth displacement fields whilst preserving the discontinuities.


### 2.1. Registration

LetΩ
⊂
R
2 be the domain of the pixel positions x
=
(
x
1
,
x
2
). We then define by the functions R
:
Ω
→
R and T
:
Ω
→
R our reference and template image. The aim of image registration is to find a transformation Φ
(
x
)
≔
x
+
w
(
x
) such that the relation T
∘
Φ
≈
R holds. The function (1)
w
:
Ω
⟶
R
2
x
⟼
w
x
≔
u
x
,
v
x
,with u
,
v
:
Ω
→
R, describes the displacement field and will be the intrinsic function we investigate. For convenience we will use the abbreviations w, u, and v for w
(
x
), u
(
x
), and v
(
x
).To solve a nonrigid registration problem with expected discontinuities in the displacement field, we want to follow a variational approach and we therefore first consider the energy functional(2)
E
f
s
w
=
∫
Ω
f
w
+
μ
s
w
d
x
,where μ
∈
R
+ is a weighting parameter and f and s are the fidelity term and the smoothness term, respectively, which are defined as(3)
f
w
=
f
u
,
v
≔
γ
1
T
x
+
w
-
R
x
+
γ
2
∂
x
1
T
x
+
w
-
∂
x
1
R
x
+
γ
2
∂
x
2
T
x
+
w
-
∂
x
2
R
x
,
(4)
s
w
=
s
u
,
v
≔
∇
u
2
+
∇
v
2
=
∇
w
.The fidelity term f incorporates the constraints for the grey value constancy and the gradient constancy and the corresponding weighting parameters are given by γ
1
,
γ
2
∈
R
0
+. The smoothness term s results in the L
1-norm and, respectively, the vectorial TV of w.The energy functionalE
f
s is motivated by the energy functional proposed by Brox et al. [8]. Instead of using the approximation function Ψ for the L
1 norm, here the pure L
1 norm is used for the smoothness term s and for the fidelity term f the sum of the absolute grey value difference and the componentwise absolute gradient differences is taken.


### 2.2. Motion Segmentation

Now we would like to integrate the formulation of the energy functional above into the convex segmentation model of Chan et al. [14]. To differentiate the displacement field w into w
+ and w
-, we therefore choose a binary function (5)
u
~
:
R
2
⟶
0,1
x
⟼
u
~
x
≔
1
Σ
x
,where Σ
⊆
Ω
⊆
R
2, with Σ
≔
x
∈
Ω
∣
u
~
x
=
1. By defining D
(
w
)
≔
f
(
w
)
+
μ
s
(
w
) as a data term, we formulate our energy functional as(6)
E
w
+
,
w
-
,
u
~
=
∫
Ω
D
w
+
u
~
x
d
x
+
∫
Ω
D
w
-
1
-
u
~
x
d
x
+
ν
∫
Ω
∇
u
~
x
d
x
,where the last term of the above energy is regularization defined by the TV norm and ν is a constant parameter.The energy functionalE (6) is very much related to the energy functional proposed by Amiaz and Kiryati [13]. Instead of applying the Heaviside function to a level set function, the binary function u
~ is used here. Furthermore the approximation function Ψ for the L
1 norm, which was used in [13], is omitted here. Instead the pure L
1 norm is used for the smoothness terms s
(
w
±
) (4) and the fidelity terms f
(
w
±
) are replaced by a sum of absolute values of the constraints (3).As pointed out by Chan et al. in [14], (6) is strongly related to the Mumford-Shah functional [4] and can be written as(7)
E
~
w
+
,
w
-
,
Σ
=
∫
Σ
D
w
+
d
x
+
∫
Ω
∖
Σ
D
w
-
d
x
+
ν
Per
⁡
Σ
;
Ω
,where Per
⁡
Σ
;
Ω denotes the perimeter of the set Σ
⊆
Ω.Remark 1.
One can show that a global minimizerΣ
m
i
n of E
~
(
w
+
,
w
-
,
·
) can be found by solving the convexified problem m
i
n
0
≤
u
~
≤
1
E
(
w
+
,
w
-
,
u
~
) and finally setting Σ
=
Σ
(
η
)
≔
x
∈
Ω
∣
u
~
(
x
)
≥
η for almost every η with η
∈
0,1 (see Proposition  1 in [16]). In fact, this holds for any data terms D
(
w
+
) and D
(
w
-
) that are measurable.Finally, we obtain the aimed displacement field by settingw
≔
w
+
u
~
+
w
-
(
1
-
u
~
).In Figure1 we illustrate our method by an example. The motion segmentation function u
~ splits the displacement field w into the two displacement fields w
+ and w
-.Figure 1
From (a) to (d): the reference and template image, the motion segmentation functionu
~, and the displacement field w.
(a)
(b)
(c)
(d)


## 2.1. Registration

LetΩ
⊂
R
2 be the domain of the pixel positions x
=
(
x
1
,
x
2
). We then define by the functions R
:
Ω
→
R and T
:
Ω
→
R our reference and template image. The aim of image registration is to find a transformation Φ
(
x
)
≔
x
+
w
(
x
) such that the relation T
∘
Φ
≈
R holds. The function (1)
w
:
Ω
⟶
R
2
x
⟼
w
x
≔
u
x
,
v
x
,with u
,
v
:
Ω
→
R, describes the displacement field and will be the intrinsic function we investigate. For convenience we will use the abbreviations w, u, and v for w
(
x
), u
(
x
), and v
(
x
).To solve a nonrigid registration problem with expected discontinuities in the displacement field, we want to follow a variational approach and we therefore first consider the energy functional(2)
E
f
s
w
=
∫
Ω
f
w
+
μ
s
w
d
x
,where μ
∈
R
+ is a weighting parameter and f and s are the fidelity term and the smoothness term, respectively, which are defined as(3)
f
w
=
f
u
,
v
≔
γ
1
T
x
+
w
-
R
x
+
γ
2
∂
x
1
T
x
+
w
-
∂
x
1
R
x
+
γ
2
∂
x
2
T
x
+
w
-
∂
x
2
R
x
,
(4)
s
w
=
s
u
,
v
≔
∇
u
2
+
∇
v
2
=
∇
w
.The fidelity term f incorporates the constraints for the grey value constancy and the gradient constancy and the corresponding weighting parameters are given by γ
1
,
γ
2
∈
R
0
+. The smoothness term s results in the L
1-norm and, respectively, the vectorial TV of w.The energy functionalE
f
s is motivated by the energy functional proposed by Brox et al. [8]. Instead of using the approximation function Ψ for the L
1 norm, here the pure L
1 norm is used for the smoothness term s and for the fidelity term f the sum of the absolute grey value difference and the componentwise absolute gradient differences is taken.


## 2.2. Motion Segmentation

Now we would like to integrate the formulation of the energy functional above into the convex segmentation model of Chan et al. [14]. To differentiate the displacement field w into w
+ and w
-, we therefore choose a binary function (5)
u
~
:
R
2
⟶
0,1
x
⟼
u
~
x
≔
1
Σ
x
,where Σ
⊆
Ω
⊆
R
2, with Σ
≔
x
∈
Ω
∣
u
~
x
=
1. By defining D
(
w
)
≔
f
(
w
)
+
μ
s
(
w
) as a data term, we formulate our energy functional as(6)
E
w
+
,
w
-
,
u
~
=
∫
Ω
D
w
+
u
~
x
d
x
+
∫
Ω
D
w
-
1
-
u
~
x
d
x
+
ν
∫
Ω
∇
u
~
x
d
x
,where the last term of the above energy is regularization defined by the TV norm and ν is a constant parameter.The energy functionalE (6) is very much related to the energy functional proposed by Amiaz and Kiryati [13]. Instead of applying the Heaviside function to a level set function, the binary function u
~ is used here. Furthermore the approximation function Ψ for the L
1 norm, which was used in [13], is omitted here. Instead the pure L
1 norm is used for the smoothness terms s
(
w
±
) (4) and the fidelity terms f
(
w
±
) are replaced by a sum of absolute values of the constraints (3).As pointed out by Chan et al. in [14], (6) is strongly related to the Mumford-Shah functional [4] and can be written as(7)
E
~
w
+
,
w
-
,
Σ
=
∫
Σ
D
w
+
d
x
+
∫
Ω
∖
Σ
D
w
-
d
x
+
ν
Per
⁡
Σ
;
Ω
,where Per
⁡
Σ
;
Ω denotes the perimeter of the set Σ
⊆
Ω.Remark 1.
One can show that a global minimizerΣ
m
i
n of E
~
(
w
+
,
w
-
,
·
) can be found by solving the convexified problem m
i
n
0
≤
u
~
≤
1
E
(
w
+
,
w
-
,
u
~
) and finally setting Σ
=
Σ
(
η
)
≔
x
∈
Ω
∣
u
~
(
x
)
≥
η for almost every η with η
∈
0,1 (see Proposition  1 in [16]). In fact, this holds for any data terms D
(
w
+
) and D
(
w
-
) that are measurable.Finally, we obtain the aimed displacement field by settingw
≔
w
+
u
~
+
w
-
(
1
-
u
~
).In Figure1 we illustrate our method by an example. The motion segmentation function u
~ splits the displacement field w into the two displacement fields w
+ and w
-.Figure 1
From (a) to (d): the reference and template image, the motion segmentation functionu
~, and the displacement field w.
(a)
(b)
(c)
(d)


## 3. Optimisation

### 3.1. Iterative Scheme

To facilitate the optimisation procedure we replace the fidelity termf in (3) by its linearised version(8)
f
w
=
γ
1
ρ
1
w
+
γ
2
ρ
2
1
w
+
γ
2
ρ
2
2
w
,where(9)
ρ
1
w
=
ρ
1
x
,
w
0
,
w
≔
T
x
+
w
0
+
∇
T
x
+
w
0
T
w
-
w
0
-
R
x
,
ρ
2
1
w
=
ρ
2
1
x
,
w
0
,
w
≔
∂
x
1
T
x
+
w
0
+
∂
x
1
x
1
T
x
+
w
0
∂
x
2
x
1
T
x
+
w
0
T
w
-
w
0
-
∂
x
1
R
x
,
ρ
2
2
w
=
ρ
2
2
x
,
w
0
,
w
≔
∂
x
2
T
x
+
w
0
+
∂
x
1
x
2
T
x
+
w
0
∂
x
2
x
2
T
x
+
w
0
T
w
-
w
0
-
∂
x
2
R
x
,with w
0 fixed.The minimization of the energy functionalE (6) with respect to w
+, w
-, and u
~ is then performed by the following iterative scheme:(1)
For fixedw
+ and w
-, solve(10)
min
u
~
∈
0,1
⁡
∫
Ω
D
w
+
u
~
x
d
x
+
∫
Ω
D
w
-
1
-
u
~
x
d
x
+
∫
Ω
ν
∇
u
~
x
d
x
.
(2)
For fixedu
~, solve(11)
min
w
+
⁡
∫
Ω
D
w
+
u
~
x
d
x
.
(3)
For fixedu
~, solve(12)
min
w
-
⁡
∫
Ω
D
w
-
1
-
u
~
x
d
x
.Although problem (10) is convex, it is important to note that the overall minimization of the energy functional E (6) is a nonconvex optimisation problem. For the minimization of nonconvex energy functionals there exist convex relaxation methods, which are able to provide solutions close to a global minimum. Recently, Strekalovskiy et al. [22] proposed such a method for nonconvex vector-valued labelling problems. In this paper however we will not make use of these kinds of methods.Note that, compared to our previous work [16], this time we did not introduce an auxiliary v
~ in the energy functional. Therefore our iterative scheme consists of only 3 steps instead of 4. The reason for this change is that we intend to use a different numerical approach to solve our problem. More precisely, to solve subproblems (10), (11), and (12) in a fast and efficient way, we follow a primal-dual approach as described by Chambolle and Pock in [17]. We therefore recapitulate in the next section the basic notations and formulations.


### 3.2. Basic Framework for the Primal-Dual Approach of Chambolle and Pock

First, we define byX and Y two finite-dimensional real vector spaces. Their inner products are denoted by ·
,
·
X and ·
,
·
Y, respectively, and their induced norms are given by ·
X
=
·
,
·
X and ·
Y
=
〈
·
,
·
〉
Y, respectively. The general nonlinear primal problem we consider is of the form(13)
min
x
∈
X
⁡
F
K
x
+
G
x
,where F
:
Y
→
0
,
+
∞ and G
:
X
→
0
,
+
∞ are proper, convex, and lower semicontinuous and the map K
:
X
→
Y is a continuous linear operator. The corresponding primal-dual formulation of (13) is the saddle-point problem(14)
min
x
∈
X
⁡
max
y
∈
Y
⁡
K
x
,
y
Y
+
G
x
-
F
∗
y
,with F
∗
:
Y
→
R
∪
+
∞ being the convex conjugate of F. We assume that the problems above have at least one solution (
x
^
,
y
^
)
∈
X
×
Y and therefore the following holds: (15)
K
x
^
∈
∂
F
∗
y
^
,
-
K
∗
y
^
∈
∂
G
x
^
,where ∂
F
∗
(
y
^
) and ∂
G
(
x
^
) are the subdifferentials of the convex functions F
∗ at y
^ and G at x
^. Furthermore we assume that F and G are “simple”; that is, the resolvent operators I
+
σ
∂
F
∗
-
1 and I
+
τ
∂
G
-
1 are easy to compute. For a convex function f the resolvent of the operator τ
∂
f at x
~ can be calculated in our case by(16)
x
=
I
+
τ
∂
f
-
1
x
~
=
arg
⁡
min
x
⁡
x
-
x
~
2
2
τ
+
f
x
.In this paper we will only make use of Algorithm  1 in [17] with the extrapolation parameter θ
=
1. Although interesting, the usage of the other proposed algorithms is left for the moment for later research.To apply Algorithm  1 in [17] to the minimization problems (10), (11), and (12), we first need to rewrite them in their discretized version, then identify the functions F and G, and finally derive the resolvent operators I
+
σ
∂
F
∗
-
1 and I
+
τ
∂
G
-
1.For the discrete setting we therefore define by(17)
x
i
,
j
=
x
1
i
,
j
,
x
2
i
,
j
=
i
h
,
j
h
,
i
=
1
,
…
,
M
,
j
=
1
,
…
,
N
,the pixel positions in the image domain with h being the spatial step size. For the calculations of the finite differences, the discrete divergence operator, the discretized inner products, and further details, we refer the reader to [17] and the references therein.In the following sections we will derive the resolvent operators for the three given minimization problems (10), (11), and (12).


### 3.3. Resolvent Operators for Problem (10)

Let us consider the continuous problem (10). As mentioned already before, we need to rewrite the problem in its discretized version to be able to apply a primal-dual approach for the minimization. We use the spaces X
=
R
M
N and Y
=
X
×
X and their inner products:(18)
s
,
t
X
=
∑
i
,
j
s
i
,
j
t
i
,
j
,
s
,
t
∈
X
,
(19)
p
,
q
Y
=
∑
i
,
j
p
i
,
j
1
q
i
,
j
1
+
p
i
,
j
2
q
i
,
j
2
,
p
=
p
1
,
p
2
,
q
=
q
1
,
q
2
∈
Y
.Using the rectangle rule we can rewrite problem (10) in the discretized version(20)
h
2
min
u
~
∈
0,1
MN
⁡
ν
∇
u
~
1
+
∑
i
,
j
D
w
i
,
j
+
u
~
i
,
j
+
D
w
i
,
j
-
1
-
u
~
i
,
j
,where the factor h
2 can be neglected, since it has no influence on the optimal solution. The norm ∇
s
1, with s
∈
X, is defined by ∇
s
1
=
∑
i
,
j
∇
s
i
,
j, where ∇
s
i
,
j
=
∇
s
i
,
j
1
2
+
∇
s
i
,
j
2
2.Comparing (20) to (13), we see that K
=
∇, F
(
∇
u
~
)
=
ν
∇
u
~
1, and G
(
u
~
)
=
∑
i
,
j
D
(
w
i
,
j
+
)
u
~
i
,
j
+
D
(
w
i
,
j
-
)
(
1
-
u
~
i
,
j
). The solution of the resolvent operator with respect to F
∗ can be derived as (21)
p
=
I
+
σ
∂
F
∗
-
1
p
~
⟹
p
i
,
j
=
ν
p
~
i
,
j
max
⁡
ν
,
p
~
i
,
j and the one with respect to G as (22)
u
~
=
I
+
τ
∂
G
-
1
u
^
⟹
u
~
i
,
j
=
min
⁡
max
⁡
u
^
i
,
j
-
τ
D
w
i
,
j
+
-
D
w
i
,
j
-
,
0
,
1
.Remark 2.
The primal-dual formulation of problem (20) is a special case of the general problem considered in Pock et al.’s work [23] and the resulting primal-dual algorithm is then the same as described there. Namely, the dual variable p is projected onto the convex set, a disc with radius ν, and with a truncation of the primal variable u
~ the projection on the feasible convex set 0,1 is achieved.Remark 3.
Observant readers probably noticed that we slightly misused the mathematical notation to make the formulas more pleasing to the eye. More precisely, instead of writingD
(
w
i
,
j
±
) we should have been using D
(
w
±
)
i
,
j, since discretization is performed after the function D is applied. In the following we will also make use of this sloppy notation for the functions ρ
1, ρ
2
(
1
), and ρ
2
(
2
) given in (9).


### 3.4. Resolvent Operators for Problem (11)

Now we consider the continuous problem (11). If we have a closer look at this minimization problem, we see that we only receive information about w
+ on the domain Σ where u
~ will be set to 1. Although theoretically reasonable, for the numerical calculations the implementation gets facilitated by having a smooth extension of w
+ to the domain Ω
∖
Σ. We therefore consider instead the problem(23)
min
w
+
⁡
∫
Ω
f
w
+
u
~
x
+
μ
s
w
+
d
x
.Comparing (11) to (23) the only difference is that the factor u
~ is not applied to the smoothness term s anymore. For the primal-dual approach we will use this time the spaces X
=
R
M
N
×
R
M
N and Y
=
X
×
X. The inner product of X is the same as in (19) and the one for Y is defined by(24)
p
,
q
Y
=
∑
i
,
j
p
i
,
j
1
q
i
,
j
1
+
p
i
,
j
2
q
i
,
j
2
+
p
i
,
j
3
q
i
,
j
3
+
p
i
,
j
4
q
i
,
j
4
,
p
=
p
1
,
p
2
,
p
3
,
p
4
,
q
=
q
1
,
q
2
,
q
3
,
q
4
∈
Y
.After the discretization of problem (23) we obtain(25)
min
w
+
∈
X
⁡
μ
∇
w
+
1
+
∑
i
,
j
γ
1
ρ
1
w
i
,
j
+
+
γ
2
ρ
2
1
w
i
,
j
+
+
γ
2
ρ
2
2
w
i
,
j
+
u
~
i
,
j
,where this time the norm ∇
s
1, with s
=
(
s
1
,
s
2
)
∈
X, is defined by ∇
s
1
=
∑
i
,
j
∇
s
i
,
jwith ∇
s
i
,
j
=
∇
s
1
i
,
j
1
2
+
∇
s
1
i
,
j
2
2
+
∇
s
2
i
,
j
1
2
+
∇
s
2
i
,
j
2
2.Comparing (25) to (13), we see that the corresponding functions are K
=
∇, F
(
∇
w
+
)
=
μ
∇
w
+
1, and(26)
G
w
+
=
∑
i
,
j
γ
1
ρ
1
w
i
,
j
+
+
γ
2
ρ
2
1
w
i
,
j
+
+
γ
2
ρ
2
2
w
i
,
j
+
u
~
i
,
j
.From the resolvent operator with respect to F
∗ we obtain(27)
q
=
I
+
σ
∂
F
∗
-
1
q
~
⟹
q
i
,
j
=
μ
q
~
i
,
j
max
⁡
μ
,
q
~
i
,
j
.The derivation of the resolvent operator with respect to G is not that straightforward and more effort has to be put in to find a suitable solution. It is common that the resolvent operators of functions, which are sums of quadratic and absolute norms, lead to so-called thresholding schemes [24, 25]. This will be also the case here. Having a closer look at the definition of G (26) and (16), we see that we have to solve(28)
w
+
=
I
+
τ
∂
G
-
1
w
0
+
=
arg
⁡
min
w
+
∈
X
⁡
w
+
-
w
0
+
X
2
2
τ
+
G
w
+
=
arg
⁡
min
w
+
=
u
+
,
v
+
∈
X
⁡
1
2
τ
∑
i
,
j
u
i
,
j
+
-
u
0
i
,
j
+
2
+
v
i
,
j
+
-
v
0
i
,
j
+
2
+
∑
i
,
j
γ
1
ρ
1
w
i
,
j
+
+
γ
2
ρ
2
1
w
i
,
j
+
+
γ
2
ρ
2
2
w
i
,
j
+
u
~
i
,
j
.For u
~
i
,
j
=
0 we can conclude from (28) that w
i
,
j
+
=
w
0
i
,
j
+. On the other hand, for u
~
i
,
j
=
1 we have to distinguish the cases(29)
ρ
1
w
i
,
j
+
>
<
=
0
,
ρ
2
1
w
i
,
j
+
>
<
=
0
,
ρ
2
2
w
i
,
j
+
>
<
=
0
,which turn out to be 27 in total. The chosen numbering for the different cases is shown in Table 1.Table 1
This table indicates the numbering of the 27 cases (see (29)).
Case number
ρ
1
(
w
i
,
j
+
)
⋄
0
ρ
2
(
1
)
(
w
i
,
j
+
)
⋄
0
ρ
2
(
2
)
(
w
i
,
j
+
)
⋄
0
1
>
>
>
2
<
>
>
3
=
>
>
4
>
<
>
5
<
<
>
6
=
<
>
7
>
>
<
8
<
>
<
9
=
>
<
10
>
<
<
11
<
<
<
12
=
<
<
13
>
=
>
14
<
=
>
15
=
=
>
16
>
=
<
17
<
=
<
18
=
=
<
19
>
=
=
20
<
=
=
21
=
=
=
22
>
>
=
23
<
>
=
24
=
>
=
25
>
<
=
26
<
<
=
27
=
<
=


#### 3.4.1. Geometric Interpretation of Problem (28) with u
~
i
,
j
=
1

Before explaining the derivation of the explicit solutions and the needed reformulation of the conditions for the different cases in (29) in more detail, we want to give a deeper insight of the geometric interpretation of problem (28) with u
~
i
,
j
=
1. To facilitate the notation in the following, we introduce the terms(30)
a
i
,
j
=
∇
T
x
+
w
0
+
i
,
j
,
b
i
,
j
=
∂
x
1
x
1
T
x
+
w
0
+
∂
x
2
x
1
T
x
+
w
0
+
i
,
j
,
c
i
,
j
=
∂
x
1
x
2
T
x
+
w
0
+
∂
x
2
x
2
T
x
+
w
0
+
i
,
j
,
(31)
x
i
,
j
=
τ
γ
1
a
i
,
j
,
y
i
,
j
=
τ
γ
2
b
i
,
j
,
z
i
,
j
=
τ
γ
2
c
i
,
j
.Using the definitions in (9) and the notation in (30), we can rewrite ρ
1
(
w
i
,
j
+
), ρ
2
(
1
)
(
w
i
,
j
+
), and ρ
2
(
2
)
(
w
i
,
j
+
) as (32)
ρ
1
w
i
,
j
+
=
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
,
ρ
2
1
w
i
,
j
+
=
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
,
ρ
2
2
w
i
,
j
+
=
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
.Note that with ρ
1
(
w
i
,
j
+
)
=
0 a line l
a
i
,
j is defined with its normal vector being parallel to the vector a
i
,
j. Similarly, ρ
2
(
1
)
(
w
i
,
j
+
)
=
0 defines a line l
b
i
,
j and ρ
2
(
2
)
(
w
i
,
j
+
)
=
0 a line l
c
i
,
j.Now we can argue similar to Zach et al. in [25] for the geometric interpretation. The mathematical structure of their minimization problem for which they derive the thresholding scheme is very similar to the one we have in (28) with u
~
i
,
j
=
1. The first term in (28), u
i
,
j
+
-
u
0
i
,
j
+
2
+
v
i
,
j
+
-
v
0
i
,
j
+
2, is the squared distance of w
i
,
j
+ to w
0
i
,
j
+, and the terms ρ
1
w
i
,
j
+, ρ
2
1
w
i
,
j
+, and ρ
2
2
w
i
,
j
+ define the unsigned distances to the lines l
a
i
,
j, l
b
i
,
j, and l
c
i
,
j, respectively. Considering now all w
δ
i
,
j
+ with a fixed distance δ to w
0
i
,
j
+, we see that problem (28) with u
~
i
,
j
=
1 is minimized for w
δ
i
,
j
+ closest to the three lines l
a
i
,
j, l
b
i
,
j, and l
c
i
,
j. See Figure 2 for an illustration.Figure 2
The geometric interpretation of minimization problem (28) with u
~
i
,
j
=
1.


#### 3.4.2. Derivation of the Explicit Solutions for Problem (28) with u
~
i
,
j
=
1

We are ready now to derive for each of the cases an explicit solution by using (28). We show the derivation of the solutions only for four different cases. The explicit solutions for the remaining cases can be calculated in a similar fashion as the presented ones.Case  1 and Similar Ones. Let us consider the first case in Table 1; that is,(33)
ρ
1
w
i
,
j
+
=
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
>
0
,
ρ
2
1
w
i
,
j
+
=
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
>
0
,
ρ
2
2
w
i
,
j
+
=
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
>
0
.Solving problem (28) for this case leads to the equation (34)
1
2
τ
2
w
i
,
j
+
-
w
0
i
,
j
+
+
γ
1
a
i
,
j
+
γ
2
b
i
,
j
+
γ
2
c
i
,
j
=
!
0
,and by using the notation in (31) the explicit solution can then be written as(35)
w
i
,
j
+
=
w
0
i
,
j
+
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
.The derivation of the explicit solutions for cases 2, 4, 5, 7, 8, 10, and 11 in Table 1 is performed similarly. There are however other cases left which need another treatment, as, for example, case number 3.Case  3 and Similar Ones. The condition for case 3 in Table 1 is given by (36)
ρ
1
w
i
,
j
+
=
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
=
0
,
ρ
2
1
w
i
,
j
+
=
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
>
0
,
ρ
2
2
w
i
,
j
+
=
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
>
0
.Here we have to distinguish the two situations a
i
,
j
≠
0 and a
i
,
j
=
0. (i)
Ifa
i
,
j
≠
0, we have the following.Fromρ
1
(
w
i
,
j
+
)
=
0 we see that the solution w
i
,
j
+ has to lie on the line l
a
i
,
j and of course the distance of it to the other lines l
b
i
,
j and l
c
i
,
j should be minimal. Therefore we can assume that the solution is of the form(37)
w
i
,
j
+
=
w
0
i
,
j
+
+
α
1
a
i
,
j
+
β
0
p
p
n
a
i
,
j
,where n
a
i
,
j is the normal to a
i
,
j. Geometrically this means that w
0
i
,
j
+ is first orthogonally projected to the line l
a
i
,
j and afterwards moved along the line such that it minimizes the distance to l
b
i
,
j and l
c
i
,
j. The unknown α
1 is determined by replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in the equation ρ
1
(
w
i
,
j
+
)
=
0 using (37), which leads to(38)
α
1
=
-
ρ
1
w
0
i
,
j
+
a
i
,
j
2
.Solving problem (28) for the current case and replacing again the term (
w
i
,
j
+
-
w
0
i
,
j
+
) using (37) lead after some simple calculations to the second parameter (39)
β
0
p
p
=
n
a
i
,
j
T
-
y
i
,
j
-
z
i
,
j
n
a
i
,
j
2
.
(ii)
Ifa
i
,
j
=
0, we have the following.Sincea
i
,
j vanishes, there is also no line l
a
i
,
j that we have to consider and we can directly derive the solution w
i
,
j
+ by just solving (28), similar to the first case in Table 1, and we get the solution (40)
w
i
,
j
+
=
w
0
i
,
j
+
-
y
i
,
j
-
z
i
,
j
.Finally, we can derive the explicit solution for cases 6, 9, 12, 13, 14, 16, 17, 22, 23, 25, and 26 in Table 1 similar to case 3 above.Case  15 and Similar Ones. Let us now consider case number 15, which was not covered so far. The condition is this time given by (41)
ρ
1
w
i
,
j
+
=
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
=
0
,
ρ
2
1
w
i
,
j
+
=
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
=
0
,
ρ
2
2
w
i
,
j
+
=
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
>
0
,and we see that we have now two equal signs, namely, for ρ
1
(
w
i
,
j
+
) and ρ
2
(
1
)
(
w
i
,
j
+
). This time we have to distinguish the situations a
i
,
j
≠
0, a
i
,
j
=
0
∧
b
i
,
j
≠
0, and a
i
,
j
=
0
∧
b
i
,
j
=
0. (i)
Ifa
i
,
j
≠
0, we have the following.Fromρ
1
(
w
i
,
j
+
)
=
0 and ρ
2
(
1
)
(
w
i
,
j
+
)
=
0 we see that the solution w
i
,
j
+ has to lie on the lines l
a
i
,
j and l
b
i
,
j and should have a minimal distance to the line l
c
i
,
j. Thus, we assume that the explicit solution is of the form(42)
w
i
,
j
+
=
w
0
i
,
j
+
+
α
1
a
i
,
j
+
β
00
p
n
a
i
,
j
.Similar to the last paragraph, replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in the equation ρ
1
(
w
i
,
j
+
)
=
0 by making use of (42) leads to the same α
1 as in (38). For the parameter β
00
p we get(43)
β
00
p
=
-
ρ
2
1
w
0
i
,
j
+
-
α
1
b
i
,
j
T
a
i
,
j
b
i
,
j
T
n
a
i
,
j
,
if
b
i
,
j
T
n
a
i
,
j
≠
0
,
n
a
i
,
j
T
-
z
i
,
j
n
a
i
,
j
2
,
if
b
i
,
j
T
n
a
i
,
j
=
0
.The case where b
i
,
j
T
n
a
i
,
j
≠
0 in (43), which is equivalent to a
i
,
j
∦
b
i
,
j, indicates that the lines l
a
i
,
j and l
b
i
,
j are not parallel and that they therefore intersect each other in one point. The solution w
i
,
j
+ we look for will be at this intersection and the corresponding β
00
p is calculated by replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in the equation ρ
2
(
1
)
(
w
i
,
j
+
)
=
0 using (42). For the case where b
i
,
j
T
n
a
i
,
j
=
0 in (43), which is equivalent to a
i
,
j
∥
b
i
,
j, we conclude that the lines l
a
i
,
j and l
b
i
,
j are parallel. This means that if the solution w
i
,
j
+ lies on the line l
a
i
,
j it necessarily also lies on the line l
b
i
,
j and the parameter β
00
p is calculated by solving problem (28) for the current case and replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) using (42).(ii)
Ifa
i
,
j
=
0
∧
b
i
,
j
≠
0, we have the following.Sincea
i
,
j vanishes, there is no need to consider the line l
a
i
,
j and we focus instead on the line l
b
i
,
j on which the solution w
i
,
j
+ has to lie too. We therefore assume the solution has the form(44)
w
i
,
j
+
=
w
0
i
,
j
+
+
α
2
b
i
,
j
+
β
00
p
n
b
i
,
j
.Replacing this time the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in the equation ρ
2
(
1
)
(
w
i
,
j
+
)
=
0 with the help of (44) results in(45)
α
2
=
-
ρ
2
1
w
0
i
,
j
+
b
i
,
j
2
.After solving problem (28) and replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) using (44) we get (46)
β
00
p
=
n
b
i
,
j
T
-
z
i
,
j
n
b
i
,
j
2
.
(iii)
Ifa
i
,
j
=
0
∧
b
i
,
j
=
0, we have the following.This time botha
i
,
j and b
i
,
j vanish and the corresponding lines l
a
i
,
j and l
b
i
,
j do not play a role in the derivation of the solution w
i
,
j
+ anymore. We can directly solve problem (28) without any prior assumption, similarly as done before in other cases, and we get (47)
w
i
,
j
+
=
w
0
i
,
j
+
-
z
i
,
j
.The derivation of the explicit solutions for cases 18, 19, 20, 24, and 27 in Table 1 is done similarly to case 15 explained above.Case  21. There is now one case left, which was not discussed so far. The corresponding condition from Table 1 is (48)
ρ
1
w
i
,
j
+
=
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
=
0
,
ρ
2
1
w
i
,
j
+
=
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
=
0
,
ρ
2
2
w
i
,
j
+
=
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
=
0
.This time we have three equal signs and we have to distinguish the situations a
i
,
j
≠
0, a
i
,
j
=
0
∧
b
i
,
j
≠
0, a
i
,
j
=
0
∧
b
i
,
j
=
0
∧
c
i
,
j
≠
0, and a
i
,
j
=
0
∧
b
i
,
j
=
0
∧
c
i
,
j
=
0. (i)
Ifa
i
,
j
≠
0, we have the following.Sinceρ
1
(
w
i
,
j
+
)
=
0, ρ
2
(
1
)
(
w
i
,
j
+
)
=
0, and ρ
2
(
2
)
(
w
i
,
j
+
)
=
0, we see that the solution w
i
,
j
+ has to lie this time on all three lines l
a
i
,
j, l
b
i
,
j, and l
c
i
,
j. For a
i
,
j
≠
0 we assume that the explicit solution is of the form(49)
w
i
,
j
+
=
w
0
i
,
j
+
+
α
1
a
i
,
j
+
β
000
n
a
i
,
j
.Like in the paragraphs before, replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in ρ
1
(
w
i
,
j
+
)
=
0 by making use of (49) leads to α
1 being equal to the one in (38). To determine β
000 in (49), we argue as follows: (a)
Ifb
i
,
j
T
n
a
i
,
j
≠
0, which is equivalent to a
i
,
j
∦
b
i
,
j, we use (49) in ρ
2
(
1
)
(
w
i
,
j
+
)
=
0 and get (50)
β
000
=
-
ρ
2
1
w
0
i
,
j
+
-
α
1
b
i
,
j
T
a
i
,
j
b
i
,
j
T
n
a
i
,
j
.
(b)
Ifb
i
,
j
T
n
a
i
,
j
=
0 and c
i
,
j
T
n
a
i
,
j
≠
0, that is, a
i
,
j
∥
b
i
,
j and a
i
,
j
∦
c
i
,
j, we use (49) in ρ
2
(
2
)
(
w
i
,
j
+
)
=
0 and get (51)
β
000
=
-
ρ
2
2
w
0
i
,
j
+
-
α
1
c
i
,
j
T
a
i
,
j
c
i
,
j
T
n
a
i
,
j
.
(c)
Ifb
i
,
j
T
n
a
i
,
j
=
0 and c
i
,
j
T
n
a
i
,
j
=
0, this means that a
i
,
j
∥
b
i
,
j and a
i
,
j
∥
c
i
,
j. Thus all the three lines l
a
i
,
j, l
b
i
,
j, and l
c
i
,
j are parallel and by solving problem (28) and using (49) we get (52)
β
000
=
0
.
(ii)
Ifa
i
,
j
=
0
∧
b
i
,
j
≠
0, we have the following.Similarly as in the discussion of (44) in case 15, we assume the solution is of the form(53)
w
i
,
j
+
=
w
0
i
,
j
+
+
α
2
b
i
,
j
+
β
000
n
b
i
,
j
.Replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in ρ
2
(
1
)
(
w
i
,
j
+
)
=
0 by making use of (53) leads to α
2 being equal to the one in (45). For β
000 in (53), we argue as follows: (a)
Ifc
i
,
j
T
n
b
i
,
j
≠
0, which is equivalent to b
i
,
j
∦
c
i
,
j, we make use of (53) in ρ
2
(
2
)
(
w
i
,
j
+
)
=
0 and get (54)
β
000
=
-
ρ
2
2
w
0
i
,
j
+
-
α
2
c
i
,
j
T
b
i
,
j
c
i
,
j
T
n
b
i
,
j
.
(b)
Ifc
i
,
j
T
n
b
i
,
j
=
0, which means b
i
,
j
∥
c
i
,
j, we get by solving (28) and using (53) (55)
β
000
=
0
.
(iii)
Ifa
i
,
j
=
0
∧
b
i
,
j
=
0
∧
c
i
,
j
≠
0, we have the following.Here,a
i
,
j and b
i
,
j vanish and we focus therefore on the line l
c
i
,
j on which the solution w
i
,
j
+ has to lie. We assume the explicit solution is of the form(56)
w
i
,
j
+
=
w
0
i
,
j
+
+
α
3
c
i
,
j
,where α
3 is calculated by using (56) to replace (
w
i
,
j
+
-
w
0
i
,
j
+
) in ρ
2
(
2
)
(
w
i
,
j
+
)
=
0 and is given by (57)
α
3
=
-
ρ
2
2
w
0
i
,
j
+
c
i
,
j
2
.
(iv)
Ifa
i
,
j
=
0
∧
b
i
,
j
=
0
∧
c
i
,
j
=
0, we have the following.In this situation we do not have to consider the linesl
a
i
,
j, l
b
i
,
j, and l
c
i
,
j at all and solving (28) results in (58)
w
i
,
j
+
=
w
0
i
,
j
+
.


#### 3.4.3. Reformulation of the Conditions for Problem (28) with u
~
i
,
j
=
1

Sincew
i
,
j
+ is unknown but w
0
i
,
j
+ in (9) is known, we have to rewrite the conditions for the 27 cases, which are defined by (29), in terms of ρ
1
(
w
0
i
,
j
+
), ρ
2
(
1
)
(
w
0
i
,
j
+
), and ρ
2
(
2
)
(
w
0
i
,
j
+
). This can be done by incorporating the derived explicit solutions and finding proper partitioning of the space spanned by ρ
1
(
w
0
i
,
j
+
), ρ
2
(
1
)
(
w
0
i
,
j
+
), and ρ
2
(
2
)
(
w
0
i
,
j
+
) (see Figure 3). A similar approach was used in the work of Zach et al. [25] to reformulate the conditions for their proposed thresholding scheme.Figure 3
The space spanned byρ
1
w
0
i
,
j
+, ρ
2
1
w
0
i
,
j
+, and ρ
2
2
w
0
i
,
j
+.As one can see in Figure3, we indicated a cube in the coordinate system that will facilitate for us the problem of finding proper partitioning of the considered space. The edge length of the cube is chosen in such a way that the critical values of the reformulated conditions are covered. The reader will understand the meaning of this sentence afterwards in a better way. For later need, we number the faces of the cube as shown in Figure 4.Figure 4
The numbering of the faces for the cube indicated in Figure3.The orientation of the face numbers in Figure4 indicates from which direction we will look at the faces.In the following we want to show how we can rewrite the conditions for the different cases. The same four groups of cases which were discussed in Section3.4.2 are considered here too and for each of them the reformulation of the conditions can be achieved in a similar fashion.Case  1 and Similar Ones. We start again with the first case in Table 1, which belongs to the group of cases that can be handled in the easiest way compared to the other ones. The condition for this case is given in (33) and by using (35) to replace the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in the condition we get (59)
ρ
1
w
i
,
j
+
=
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
,
ρ
2
1
w
i
,
j
+
=
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
,
ρ
2
2
w
i
,
j
+
=
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
,or rewritten(60)
ρ
1
w
0
i
,
j
+
>
-
a
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
ρ
2
1
w
0
i
,
j
+
>
-
b
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
ρ
2
2
w
0
i
,
j
+
>
-
c
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
.One of the critical values we mentioned before concerning the edge length of the cube shown in Figure 3 can be determined now for case number 1, namely, by the values on the right-hand side in (60).The reformulation of the conditions for cases2, 4, 5, 7, 8, 10, and 11 in Table 1 is achieved similarly to this case and the remaining seven critical values, which are necessary to define the edge length of the cube, are determined in the same way as above.Remark 4.
Before jumping to the next group of cases, we want to explain why the cube in Figure3 will help us in finding proper partitioning of the space spanned by ρ
1
(
w
0
i
,
j
+
), ρ
2
(
1
)
(
w
0
i
,
j
+
), and ρ
2
(
2
)
(
w
0
i
,
j
+
). Our aim is to divide the space into 27 parts, each one defining a region for which exactly one of the 27 cases may apply with their already derived explicit solutions. This partitioning should not be done just arbitrarily and has to be meaningful and to conform to the already set up formulations. Before, we mentioned once that we want the eight critical values from the paragraph before to be covered by the cube. The reason for this is that like this it is possible to reduce the partitioning problem of the whole space to a partitioning problem of the cube. The uniquely defined regions of the cube are the ones from the corresponding eight cases discussed in the paragraph above and they define the corners of it. In Figure 5 we show the location of the regions for different case numbers. For simplicity we drew all the regions with the same size, although they can and most likely will have different sizes. Note that the region for case number 21 is not visible in the figure, because it fully lies in the interior of the cube.
Let us focus now on one face of the cube. We know the values for the boundaries of the corner regions from the critical values we got from the previous paragraph. They are fixed and we introduce for them the namesA
v, A
h, B
v, B
h, C
v, C
h, D
v, and D
h as shown in Figure 6.
For face number6, for example, the values D
v and D
h will be the boundary values to region 1 (see Figure 5) and after having a look at (60) we can conclude that D
v
=
-
a
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j and D
h
=
-
b
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j.
For each face it always holds that(61)
A
v
≤
B
v
,
C
v
≤
D
v
,
B
h
≤
D
h
,
A
h
≤
C
h
.For example, to see that C
v
≤
D
v we consider the reformulation of the condition for case number 2:(62)
ρ
1
w
0
i
,
j
+
<
-
a
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
ρ
2
1
w
0
i
,
j
+
>
-
b
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
ρ
2
2
w
0
i
,
j
+
>
-
c
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,which can be derived, as mentioned before, similarly to case number 1. For face number 6 we see in the same manner as before by having a look at Figure 5 and (62) that C
v
=
-
a
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j. Since a
i
,
j
T
x
i
,
j
=
τ
γ
1
a
i
,
j
2
≥
0 (see (31)), we can conclude that C
v
≤
D
v.
Let us consider again an arbitrary face of the cube. Although we have boundaries for the corner regions, the ones for the regions in between are not fully defined. Since a rectangular division of the face seems to be adequate, we investigate the problem on how to achieve such one. In Figure7 a face is depicted with its corner regions and two ways of possible divisions into 9 rectangles are shown, which we call the “clockwise” and “anticlockwise” division method.
Depending on the values forA
v, A
h, B
v, B
h, C
v, C
h, D
v, and D
h sometimes both the clockwise and the anticlockwise division method work, as, for example, in Figure 7. Depending on the situation however, there can be also values for which only the clockwise or the anticlockwise or even none of both methods will work, where in the latter situation another adequate division method has to be defined.
In Table2 we list the different possible situations for the values A
v, A
h, B
v, B
h, C
v, C
h, D
v, and D
h.
Furthermore, we split the last situation S9 into its parts for later use and introduce the corresponding labelling in Table3.
In the following, we go through the different situations given in Tables2 and 3 and decide which kind of method can be used to achieve a division of the face into 9 rectangular regions: (i)
Situation S1 in Table2 is illustrated in Figure 8 for face number 6 and the corresponding region numbers. For this situation, neither the clockwise nor the anticlockwise division method will work, since the corner regions 1 and 5 overlap. An adapted division method, which is depicted in Figure 8, will be used instead for this situation. Note that in this situation certain corner regions, which were defined through the reformulation of the conditions for the group of cases similar to case 1, have to be adapted by removing from them the region where they overlap.
(ii)
For the situations S2, S3, S5, and S6 in Table2 and S9.1, S9.2, S9.3, S9.4, and S9.5 in Table 3 we can use the clockwise division method illustrated in Figure 7. We define the set of situations for which the clockwise division method can be applied by (63)
S
c
≔
S
2
,
S
3
,
S
5
,
S
6
,
S
9.1
,
S
9.2
,
S
9.3
,
S
9.4
,
S
9.5
.
(iii)
Situation S9.7 in Table3 needs again a special treatment to achieve a meaningful division, since neither the clockwise nor the anticlockwise division method can be applied because of the overlapping of the corner regions 2 and 4, if we consider face number 6, for example. The adapted division method, which will be used instead, is shown in Figure 9 for face number 6 and the corresponding region numbers. Note again that also in this situation certain corner regions have to be adapted by removing from them the region where they overlap.
(iv)
Finally, for the situations S4, S7, and S8 in Table2 and S9.6, S9.8, and S9.9 in Table 3 we can use this time the anticlockwise division method illustrated in Figure 7. The set of situations for which the anticlockwise division method can be applied is then given by(64)
S
a
c
≔
S
4
,
S
7
,
S
8
,
S
9.6
,
S
9.8
,
S
9.9
.We are now able to give an idea on how to reformulate the conditions for the remaining group of cases.Table 2
Different possible situations for the valuesA
v, A
h, B
v, B
h, C
v, C
h, D
v, and D
h.
Number
A
v
⋄
D
v
B
v
⋄
C
v
A
h
⋄
D
h
B
h
⋄
C
h
S1
>
>
>
<
S2
>
>
=
<
=
S3
>
>
<
>
<
=
S4
=
>
=
>
<
S5
=
>
=
=
<
=
S6
=
>
=
<
>
<
=
S7
<
>
<
=
>
<
S8
<
>
<
=
=
<
=
S9
<
>
<
=
<
>
<
=Table 3
Splitting of situation S9 in Table2 into its parts.
Number
A
v
⋄
D
v
B
v
⋄
C
v
A
h
⋄
D
h
B
h
⋄
C
h
S9.1
<
>
<
>
S9.2
<
>
<
=
S9.3
<
>
<
<
S9.4
<
=
<
>
S9.5
<
=
<
=
S9.6
<
=
<
<
S9.7
<
<
<
>
S9.8
<
<
<
=
S9.9
<
<
<
<Figure 5
The net of the cube indicated in Figure3 with the corresponding case numbers for the different regions. The face numbers from Figure 4 are visible in the background.Figure 6
One face of the cube with the depicted corner regions and their labelling for the boundary values. The arrows indicate the direction of increasement.Figure 7
A face with its depicted corner regions (a) and the clockwise (b) and anticlockwise (c) division method, respectively.
(a)
(b)
(c)Figure 8
The division method used for situation S1 in Table2. For the illustration we used face number 6 and the corresponding region numbers.Figure 9
The division method used for situation S9.7 in Table3. For the illustration we used again face number 6 and the corresponding region numbers.Case  3 and Similar Ones. The reformulation of the cases in this group with the same approach as for case number 1 will not lead to proper partitioning of the space shown in Figure 3. Instead we will use a more geometrically based approach which depends on the different situations the boundary values can take.First we can derive for sure, independently of which approach we use, that(65)
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
-
y
i
,
j
-
z
i
,
j
≤
a
i
,
j
T
x
i
,
j
.This can be seen by considering Figure 5, from which we can deduce that region 3 has to lie between regions 1 and 2. By using the critical values of these adjacent regions, which can be determined from (60) and (62), we arrive at the inequality above.For the reformulation of the condition with respect toρ
2
(
1
)
(
w
0
i
,
j
+
), we consider the different situations in which the boundary values of the corner regions can be for face number 6.(i)
If the boundary values are in situation S1 given in Table2, the reformulation is given by (66)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
+
x
i
,
j
+
y
i
,
j
-
z
i
,
j
≥
0
.
(ii)
If the boundary values are in one of the situations in the setS
c defined in (63) or in situation S9.7 in Table 3, we get (67)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
.
(iii)
If the boundary values are in one of the situations inS
a
c (64), the reformulation is (68)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
.For the reformulation of the condition with respect to ρ
2
(
2
)
(
w
0
i
,
j
+
), we consider this time the different situations belonging to face number 4.(i)
If the boundary values are in situation S1, the reformulation is given by(69)
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
-
x
i
,
j
-
y
i
,
j
+
z
i
,
j
≥
0
.
(ii)
If the situation of the boundary values is in the setS
c or corresponds to the situation S9.7, we get (70)
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
.
(iii)
Finally, if the situation of the boundary values is inS
a
c, the reformulation is (71)
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
.Finally, the reformulation of the conditions for cases 6, 9, 12, 13, 14, 16, 17, 22, 23, 25, and 26 in Table 1 is achieved similarly to this case.Case  15 and Similar Ones. For this group of cases we will again make use of the different situations of the boundary values, but this time only for face number 6.For the reformulation of the condition with respect toρ
1
(
w
0
i
,
j
+
), we get the following list:(i)
If the boundary values are in situation S1(72)
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
-
z
i
,
j
<
-
a
i
,
j
T
+
x
i
,
j
+
y
i
,
j
.
(ii)
If the situation of the boundary values is in the setS
c
(73)
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
-
z
i
,
j
≤
-
a
i
,
j
T
-
x
i
,
j
+
y
i
,
j
.
(iii)
If the boundary values are in situation S9.7(74)
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
-
z
i
,
j
<
-
a
i
,
j
T
+
x
i
,
j
-
y
i
,
j
.
(iv)
If the situation of the boundary values is in the setS
a
c
(75)
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
-
z
i
,
j
≤
-
a
i
,
j
T
-
x
i
,
j
-
y
i
,
j
.For the reformulation of the condition with respect toρ
2
(
1
)
(
w
0
i
,
j
+
), we get the following list:(i)
If the boundary values are in situation S1(76)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
-
z
i
,
j
<
-
b
i
,
j
T
+
x
i
,
j
+
y
i
,
j
.
(ii)
If the situation of the boundary values is in the setS
c
(77)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
-
z
i
,
j
≤
-
b
i
,
j
T
-
x
i
,
j
-
y
i
,
j
.
(iii)
If the boundary values are in situation S9.7(78)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
-
z
i
,
j
<
-
b
i
,
j
T
-
x
i
,
j
+
y
i
,
j
.
(iv)
If the situation of the boundary values is in the setS
a
c
(79)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
-
z
i
,
j
≤
-
b
i
,
j
T
+
x
i
,
j
-
y
i
,
j
.Finally, we get for the reformulation of the condition with respect to ρ
2
(
2
)
(
w
0
i
,
j
+
) the following list:(i)
If the boundary values are in situation S1(80)
ρ
2
2
w
0
i
,
j
+
+
min
⁡
c
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
c
i
,
j
T
+
x
i
,
j
+
y
i
,
j
-
z
i
,
j
>
0
.
(ii)
If the situation of the boundary values is in the setS
c or S
a
c
(81)
ρ
2
2
w
0
i
,
j
+
+
min
⁡
c
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
c
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
c
i
,
j
T
-
x
i
,
j
+
y
i
,
j
-
z
i
,
j
,
c
i
,
j
T
+
x
i
,
j
+
y
i
,
j
-
z
i
,
j
>
0
.
(iii)
If the boundary values are in situation S9.7(82)
ρ
2
2
w
0
i
,
j
+
+
min
⁡
c
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
c
i
,
j
T
-
x
i
,
j
+
y
i
,
j
-
z
i
,
j
>
0
.The reformulation of the conditions for cases 18, 19, 20, 24, and 27 in Table 1 is done then similarly to case 15 explained above.Case  21. The last case which is left is case number 21. The corresponding region contains the part of the cube which was not covered so far with the corresponding regions of the other cases. Additionally, since it is unfortunately also possible that there are overlappings for some of the regions which appear at different faces, as, for example, the possible overlapping of region 1 with region 11, these overlappings are also assigned to case number 21 to finally achieve proper partitioning.One could get the impression that this handling of the overlappings is mathematically rather grubby, but to legitimate this we recall that an optimal solution of problem (28) will result in rather small values for ρ
1
(
w
i
,
j
+
), ρ
2
(
1
)
(
w
i
,
j
+
), and ρ
2
(
2
)
(
w
i
,
j
+
) and therefore should be ideally close to case 21 in Table 1. Extending therefore the region for case number 21 should not be harmful. Furthermore, for small values of ρ
1
(
w
i
,
j
+
), ρ
2
(
1
)
(
w
i
,
j
+
), and ρ
2
(
2
)
(
w
i
,
j
+
) a reformulation of the condition can lead more likely to overlappings of certain region parts. Therefore it seems to be adequate to assign these overlappings to region 21.It is possible to define region21 in mathematical terms, but since we explained the content of this region already above this will be rather uninspiring and for the implementation an explicit formulation of the region is also not needed, since it can easily be depicted with the help of the other regions.


### 3.5. Resolvent Operators for Problem (12)

This section is very similar to Section3.4 and we mainly have to just replace w
+ by w
- and u
~ by 1
-
u
~. To achieve a smooth extension of w
- to the domain Σ, we consider now instead of problem (12) the following problem:(83)
min
w
-
⁡
∫
Ω
f
w
-
1
-
u
~
x
+
μ
s
w
-
d
x
.The spaces X and Y and their inner products are defined in the same way as in Section 3.4 and the discretization of (83) is performed in the same manner as before.The resolvent operator with respect toF
∗ is identical to (27). Finally, compared to the section before, there are only very slight changes of the resolvent operator with respect to G.


## 3.1. Iterative Scheme

To facilitate the optimisation procedure we replace the fidelity termf in (3) by its linearised version(8)
f
w
=
γ
1
ρ
1
w
+
γ
2
ρ
2
1
w
+
γ
2
ρ
2
2
w
,where(9)
ρ
1
w
=
ρ
1
x
,
w
0
,
w
≔
T
x
+
w
0
+
∇
T
x
+
w
0
T
w
-
w
0
-
R
x
,
ρ
2
1
w
=
ρ
2
1
x
,
w
0
,
w
≔
∂
x
1
T
x
+
w
0
+
∂
x
1
x
1
T
x
+
w
0
∂
x
2
x
1
T
x
+
w
0
T
w
-
w
0
-
∂
x
1
R
x
,
ρ
2
2
w
=
ρ
2
2
x
,
w
0
,
w
≔
∂
x
2
T
x
+
w
0
+
∂
x
1
x
2
T
x
+
w
0
∂
x
2
x
2
T
x
+
w
0
T
w
-
w
0
-
∂
x
2
R
x
,with w
0 fixed.The minimization of the energy functionalE (6) with respect to w
+, w
-, and u
~ is then performed by the following iterative scheme:(1)
For fixedw
+ and w
-, solve(10)
min
u
~
∈
0,1
⁡
∫
Ω
D
w
+
u
~
x
d
x
+
∫
Ω
D
w
-
1
-
u
~
x
d
x
+
∫
Ω
ν
∇
u
~
x
d
x
.
(2)
For fixedu
~, solve(11)
min
w
+
⁡
∫
Ω
D
w
+
u
~
x
d
x
.
(3)
For fixedu
~, solve(12)
min
w
-
⁡
∫
Ω
D
w
-
1
-
u
~
x
d
x
.Although problem (10) is convex, it is important to note that the overall minimization of the energy functional E (6) is a nonconvex optimisation problem. For the minimization of nonconvex energy functionals there exist convex relaxation methods, which are able to provide solutions close to a global minimum. Recently, Strekalovskiy et al. [22] proposed such a method for nonconvex vector-valued labelling problems. In this paper however we will not make use of these kinds of methods.Note that, compared to our previous work [16], this time we did not introduce an auxiliary v
~ in the energy functional. Therefore our iterative scheme consists of only 3 steps instead of 4. The reason for this change is that we intend to use a different numerical approach to solve our problem. More precisely, to solve subproblems (10), (11), and (12) in a fast and efficient way, we follow a primal-dual approach as described by Chambolle and Pock in [17]. We therefore recapitulate in the next section the basic notations and formulations.


## 3.2. Basic Framework for the Primal-Dual Approach of Chambolle and Pock

First, we define byX and Y two finite-dimensional real vector spaces. Their inner products are denoted by ·
,
·
X and ·
,
·
Y, respectively, and their induced norms are given by ·
X
=
·
,
·
X and ·
Y
=
〈
·
,
·
〉
Y, respectively. The general nonlinear primal problem we consider is of the form(13)
min
x
∈
X
⁡
F
K
x
+
G
x
,where F
:
Y
→
0
,
+
∞ and G
:
X
→
0
,
+
∞ are proper, convex, and lower semicontinuous and the map K
:
X
→
Y is a continuous linear operator. The corresponding primal-dual formulation of (13) is the saddle-point problem(14)
min
x
∈
X
⁡
max
y
∈
Y
⁡
K
x
,
y
Y
+
G
x
-
F
∗
y
,with F
∗
:
Y
→
R
∪
+
∞ being the convex conjugate of F. We assume that the problems above have at least one solution (
x
^
,
y
^
)
∈
X
×
Y and therefore the following holds: (15)
K
x
^
∈
∂
F
∗
y
^
,
-
K
∗
y
^
∈
∂
G
x
^
,where ∂
F
∗
(
y
^
) and ∂
G
(
x
^
) are the subdifferentials of the convex functions F
∗ at y
^ and G at x
^. Furthermore we assume that F and G are “simple”; that is, the resolvent operators I
+
σ
∂
F
∗
-
1 and I
+
τ
∂
G
-
1 are easy to compute. For a convex function f the resolvent of the operator τ
∂
f at x
~ can be calculated in our case by(16)
x
=
I
+
τ
∂
f
-
1
x
~
=
arg
⁡
min
x
⁡
x
-
x
~
2
2
τ
+
f
x
.In this paper we will only make use of Algorithm  1 in [17] with the extrapolation parameter θ
=
1. Although interesting, the usage of the other proposed algorithms is left for the moment for later research.To apply Algorithm  1 in [17] to the minimization problems (10), (11), and (12), we first need to rewrite them in their discretized version, then identify the functions F and G, and finally derive the resolvent operators I
+
σ
∂
F
∗
-
1 and I
+
τ
∂
G
-
1.For the discrete setting we therefore define by(17)
x
i
,
j
=
x
1
i
,
j
,
x
2
i
,
j
=
i
h
,
j
h
,
i
=
1
,
…
,
M
,
j
=
1
,
…
,
N
,the pixel positions in the image domain with h being the spatial step size. For the calculations of the finite differences, the discrete divergence operator, the discretized inner products, and further details, we refer the reader to [17] and the references therein.In the following sections we will derive the resolvent operators for the three given minimization problems (10), (11), and (12).


## 3.3. Resolvent Operators for Problem (10)

Let us consider the continuous problem (10). As mentioned already before, we need to rewrite the problem in its discretized version to be able to apply a primal-dual approach for the minimization. We use the spaces X
=
R
M
N and Y
=
X
×
X and their inner products:(18)
s
,
t
X
=
∑
i
,
j
s
i
,
j
t
i
,
j
,
s
,
t
∈
X
,
(19)
p
,
q
Y
=
∑
i
,
j
p
i
,
j
1
q
i
,
j
1
+
p
i
,
j
2
q
i
,
j
2
,
p
=
p
1
,
p
2
,
q
=
q
1
,
q
2
∈
Y
.Using the rectangle rule we can rewrite problem (10) in the discretized version(20)
h
2
min
u
~
∈
0,1
MN
⁡
ν
∇
u
~
1
+
∑
i
,
j
D
w
i
,
j
+
u
~
i
,
j
+
D
w
i
,
j
-
1
-
u
~
i
,
j
,where the factor h
2 can be neglected, since it has no influence on the optimal solution. The norm ∇
s
1, with s
∈
X, is defined by ∇
s
1
=
∑
i
,
j
∇
s
i
,
j, where ∇
s
i
,
j
=
∇
s
i
,
j
1
2
+
∇
s
i
,
j
2
2.Comparing (20) to (13), we see that K
=
∇, F
(
∇
u
~
)
=
ν
∇
u
~
1, and G
(
u
~
)
=
∑
i
,
j
D
(
w
i
,
j
+
)
u
~
i
,
j
+
D
(
w
i
,
j
-
)
(
1
-
u
~
i
,
j
). The solution of the resolvent operator with respect to F
∗ can be derived as (21)
p
=
I
+
σ
∂
F
∗
-
1
p
~
⟹
p
i
,
j
=
ν
p
~
i
,
j
max
⁡
ν
,
p
~
i
,
j and the one with respect to G as (22)
u
~
=
I
+
τ
∂
G
-
1
u
^
⟹
u
~
i
,
j
=
min
⁡
max
⁡
u
^
i
,
j
-
τ
D
w
i
,
j
+
-
D
w
i
,
j
-
,
0
,
1
.Remark 2.
The primal-dual formulation of problem (20) is a special case of the general problem considered in Pock et al.’s work [23] and the resulting primal-dual algorithm is then the same as described there. Namely, the dual variable p is projected onto the convex set, a disc with radius ν, and with a truncation of the primal variable u
~ the projection on the feasible convex set 0,1 is achieved.Remark 3.
Observant readers probably noticed that we slightly misused the mathematical notation to make the formulas more pleasing to the eye. More precisely, instead of writingD
(
w
i
,
j
±
) we should have been using D
(
w
±
)
i
,
j, since discretization is performed after the function D is applied. In the following we will also make use of this sloppy notation for the functions ρ
1, ρ
2
(
1
), and ρ
2
(
2
) given in (9).


## 3.4. Resolvent Operators for Problem (11)

Now we consider the continuous problem (11). If we have a closer look at this minimization problem, we see that we only receive information about w
+ on the domain Σ where u
~ will be set to 1. Although theoretically reasonable, for the numerical calculations the implementation gets facilitated by having a smooth extension of w
+ to the domain Ω
∖
Σ. We therefore consider instead the problem(23)
min
w
+
⁡
∫
Ω
f
w
+
u
~
x
+
μ
s
w
+
d
x
.Comparing (11) to (23) the only difference is that the factor u
~ is not applied to the smoothness term s anymore. For the primal-dual approach we will use this time the spaces X
=
R
M
N
×
R
M
N and Y
=
X
×
X. The inner product of X is the same as in (19) and the one for Y is defined by(24)
p
,
q
Y
=
∑
i
,
j
p
i
,
j
1
q
i
,
j
1
+
p
i
,
j
2
q
i
,
j
2
+
p
i
,
j
3
q
i
,
j
3
+
p
i
,
j
4
q
i
,
j
4
,
p
=
p
1
,
p
2
,
p
3
,
p
4
,
q
=
q
1
,
q
2
,
q
3
,
q
4
∈
Y
.After the discretization of problem (23) we obtain(25)
min
w
+
∈
X
⁡
μ
∇
w
+
1
+
∑
i
,
j
γ
1
ρ
1
w
i
,
j
+
+
γ
2
ρ
2
1
w
i
,
j
+
+
γ
2
ρ
2
2
w
i
,
j
+
u
~
i
,
j
,where this time the norm ∇
s
1, with s
=
(
s
1
,
s
2
)
∈
X, is defined by ∇
s
1
=
∑
i
,
j
∇
s
i
,
jwith ∇
s
i
,
j
=
∇
s
1
i
,
j
1
2
+
∇
s
1
i
,
j
2
2
+
∇
s
2
i
,
j
1
2
+
∇
s
2
i
,
j
2
2.Comparing (25) to (13), we see that the corresponding functions are K
=
∇, F
(
∇
w
+
)
=
μ
∇
w
+
1, and(26)
G
w
+
=
∑
i
,
j
γ
1
ρ
1
w
i
,
j
+
+
γ
2
ρ
2
1
w
i
,
j
+
+
γ
2
ρ
2
2
w
i
,
j
+
u
~
i
,
j
.From the resolvent operator with respect to F
∗ we obtain(27)
q
=
I
+
σ
∂
F
∗
-
1
q
~
⟹
q
i
,
j
=
μ
q
~
i
,
j
max
⁡
μ
,
q
~
i
,
j
.The derivation of the resolvent operator with respect to G is not that straightforward and more effort has to be put in to find a suitable solution. It is common that the resolvent operators of functions, which are sums of quadratic and absolute norms, lead to so-called thresholding schemes [24, 25]. This will be also the case here. Having a closer look at the definition of G (26) and (16), we see that we have to solve(28)
w
+
=
I
+
τ
∂
G
-
1
w
0
+
=
arg
⁡
min
w
+
∈
X
⁡
w
+
-
w
0
+
X
2
2
τ
+
G
w
+
=
arg
⁡
min
w
+
=
u
+
,
v
+
∈
X
⁡
1
2
τ
∑
i
,
j
u
i
,
j
+
-
u
0
i
,
j
+
2
+
v
i
,
j
+
-
v
0
i
,
j
+
2
+
∑
i
,
j
γ
1
ρ
1
w
i
,
j
+
+
γ
2
ρ
2
1
w
i
,
j
+
+
γ
2
ρ
2
2
w
i
,
j
+
u
~
i
,
j
.For u
~
i
,
j
=
0 we can conclude from (28) that w
i
,
j
+
=
w
0
i
,
j
+. On the other hand, for u
~
i
,
j
=
1 we have to distinguish the cases(29)
ρ
1
w
i
,
j
+
>
<
=
0
,
ρ
2
1
w
i
,
j
+
>
<
=
0
,
ρ
2
2
w
i
,
j
+
>
<
=
0
,which turn out to be 27 in total. The chosen numbering for the different cases is shown in Table 1.Table 1
This table indicates the numbering of the 27 cases (see (29)).
Case number
ρ
1
(
w
i
,
j
+
)
⋄
0
ρ
2
(
1
)
(
w
i
,
j
+
)
⋄
0
ρ
2
(
2
)
(
w
i
,
j
+
)
⋄
0
1
>
>
>
2
<
>
>
3
=
>
>
4
>
<
>
5
<
<
>
6
=
<
>
7
>
>
<
8
<
>
<
9
=
>
<
10
>
<
<
11
<
<
<
12
=
<
<
13
>
=
>
14
<
=
>
15
=
=
>
16
>
=
<
17
<
=
<
18
=
=
<
19
>
=
=
20
<
=
=
21
=
=
=
22
>
>
=
23
<
>
=
24
=
>
=
25
>
<
=
26
<
<
=
27
=
<
=


### 3.4.1. Geometric Interpretation of Problem (28) with u
~
i
,
j
=
1

Before explaining the derivation of the explicit solutions and the needed reformulation of the conditions for the different cases in (29) in more detail, we want to give a deeper insight of the geometric interpretation of problem (28) with u
~
i
,
j
=
1. To facilitate the notation in the following, we introduce the terms(30)
a
i
,
j
=
∇
T
x
+
w
0
+
i
,
j
,
b
i
,
j
=
∂
x
1
x
1
T
x
+
w
0
+
∂
x
2
x
1
T
x
+
w
0
+
i
,
j
,
c
i
,
j
=
∂
x
1
x
2
T
x
+
w
0
+
∂
x
2
x
2
T
x
+
w
0
+
i
,
j
,
(31)
x
i
,
j
=
τ
γ
1
a
i
,
j
,
y
i
,
j
=
τ
γ
2
b
i
,
j
,
z
i
,
j
=
τ
γ
2
c
i
,
j
.Using the definitions in (9) and the notation in (30), we can rewrite ρ
1
(
w
i
,
j
+
), ρ
2
(
1
)
(
w
i
,
j
+
), and ρ
2
(
2
)
(
w
i
,
j
+
) as (32)
ρ
1
w
i
,
j
+
=
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
,
ρ
2
1
w
i
,
j
+
=
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
,
ρ
2
2
w
i
,
j
+
=
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
.Note that with ρ
1
(
w
i
,
j
+
)
=
0 a line l
a
i
,
j is defined with its normal vector being parallel to the vector a
i
,
j. Similarly, ρ
2
(
1
)
(
w
i
,
j
+
)
=
0 defines a line l
b
i
,
j and ρ
2
(
2
)
(
w
i
,
j
+
)
=
0 a line l
c
i
,
j.Now we can argue similar to Zach et al. in [25] for the geometric interpretation. The mathematical structure of their minimization problem for which they derive the thresholding scheme is very similar to the one we have in (28) with u
~
i
,
j
=
1. The first term in (28), u
i
,
j
+
-
u
0
i
,
j
+
2
+
v
i
,
j
+
-
v
0
i
,
j
+
2, is the squared distance of w
i
,
j
+ to w
0
i
,
j
+, and the terms ρ
1
w
i
,
j
+, ρ
2
1
w
i
,
j
+, and ρ
2
2
w
i
,
j
+ define the unsigned distances to the lines l
a
i
,
j, l
b
i
,
j, and l
c
i
,
j, respectively. Considering now all w
δ
i
,
j
+ with a fixed distance δ to w
0
i
,
j
+, we see that problem (28) with u
~
i
,
j
=
1 is minimized for w
δ
i
,
j
+ closest to the three lines l
a
i
,
j, l
b
i
,
j, and l
c
i
,
j. See Figure 2 for an illustration.Figure 2
The geometric interpretation of minimization problem (28) with u
~
i
,
j
=
1.


### 3.4.2. Derivation of the Explicit Solutions for Problem (28) with u
~
i
,
j
=
1

We are ready now to derive for each of the cases an explicit solution by using (28). We show the derivation of the solutions only for four different cases. The explicit solutions for the remaining cases can be calculated in a similar fashion as the presented ones.Case  1 and Similar Ones. Let us consider the first case in Table 1; that is,(33)
ρ
1
w
i
,
j
+
=
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
>
0
,
ρ
2
1
w
i
,
j
+
=
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
>
0
,
ρ
2
2
w
i
,
j
+
=
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
>
0
.Solving problem (28) for this case leads to the equation (34)
1
2
τ
2
w
i
,
j
+
-
w
0
i
,
j
+
+
γ
1
a
i
,
j
+
γ
2
b
i
,
j
+
γ
2
c
i
,
j
=
!
0
,and by using the notation in (31) the explicit solution can then be written as(35)
w
i
,
j
+
=
w
0
i
,
j
+
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
.The derivation of the explicit solutions for cases 2, 4, 5, 7, 8, 10, and 11 in Table 1 is performed similarly. There are however other cases left which need another treatment, as, for example, case number 3.Case  3 and Similar Ones. The condition for case 3 in Table 1 is given by (36)
ρ
1
w
i
,
j
+
=
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
=
0
,
ρ
2
1
w
i
,
j
+
=
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
>
0
,
ρ
2
2
w
i
,
j
+
=
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
>
0
.Here we have to distinguish the two situations a
i
,
j
≠
0 and a
i
,
j
=
0. (i)
Ifa
i
,
j
≠
0, we have the following.Fromρ
1
(
w
i
,
j
+
)
=
0 we see that the solution w
i
,
j
+ has to lie on the line l
a
i
,
j and of course the distance of it to the other lines l
b
i
,
j and l
c
i
,
j should be minimal. Therefore we can assume that the solution is of the form(37)
w
i
,
j
+
=
w
0
i
,
j
+
+
α
1
a
i
,
j
+
β
0
p
p
n
a
i
,
j
,where n
a
i
,
j is the normal to a
i
,
j. Geometrically this means that w
0
i
,
j
+ is first orthogonally projected to the line l
a
i
,
j and afterwards moved along the line such that it minimizes the distance to l
b
i
,
j and l
c
i
,
j. The unknown α
1 is determined by replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in the equation ρ
1
(
w
i
,
j
+
)
=
0 using (37), which leads to(38)
α
1
=
-
ρ
1
w
0
i
,
j
+
a
i
,
j
2
.Solving problem (28) for the current case and replacing again the term (
w
i
,
j
+
-
w
0
i
,
j
+
) using (37) lead after some simple calculations to the second parameter (39)
β
0
p
p
=
n
a
i
,
j
T
-
y
i
,
j
-
z
i
,
j
n
a
i
,
j
2
.
(ii)
Ifa
i
,
j
=
0, we have the following.Sincea
i
,
j vanishes, there is also no line l
a
i
,
j that we have to consider and we can directly derive the solution w
i
,
j
+ by just solving (28), similar to the first case in Table 1, and we get the solution (40)
w
i
,
j
+
=
w
0
i
,
j
+
-
y
i
,
j
-
z
i
,
j
.Finally, we can derive the explicit solution for cases 6, 9, 12, 13, 14, 16, 17, 22, 23, 25, and 26 in Table 1 similar to case 3 above.Case  15 and Similar Ones. Let us now consider case number 15, which was not covered so far. The condition is this time given by (41)
ρ
1
w
i
,
j
+
=
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
=
0
,
ρ
2
1
w
i
,
j
+
=
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
=
0
,
ρ
2
2
w
i
,
j
+
=
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
>
0
,and we see that we have now two equal signs, namely, for ρ
1
(
w
i
,
j
+
) and ρ
2
(
1
)
(
w
i
,
j
+
). This time we have to distinguish the situations a
i
,
j
≠
0, a
i
,
j
=
0
∧
b
i
,
j
≠
0, and a
i
,
j
=
0
∧
b
i
,
j
=
0. (i)
Ifa
i
,
j
≠
0, we have the following.Fromρ
1
(
w
i
,
j
+
)
=
0 and ρ
2
(
1
)
(
w
i
,
j
+
)
=
0 we see that the solution w
i
,
j
+ has to lie on the lines l
a
i
,
j and l
b
i
,
j and should have a minimal distance to the line l
c
i
,
j. Thus, we assume that the explicit solution is of the form(42)
w
i
,
j
+
=
w
0
i
,
j
+
+
α
1
a
i
,
j
+
β
00
p
n
a
i
,
j
.Similar to the last paragraph, replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in the equation ρ
1
(
w
i
,
j
+
)
=
0 by making use of (42) leads to the same α
1 as in (38). For the parameter β
00
p we get(43)
β
00
p
=
-
ρ
2
1
w
0
i
,
j
+
-
α
1
b
i
,
j
T
a
i
,
j
b
i
,
j
T
n
a
i
,
j
,
if
b
i
,
j
T
n
a
i
,
j
≠
0
,
n
a
i
,
j
T
-
z
i
,
j
n
a
i
,
j
2
,
if
b
i
,
j
T
n
a
i
,
j
=
0
.The case where b
i
,
j
T
n
a
i
,
j
≠
0 in (43), which is equivalent to a
i
,
j
∦
b
i
,
j, indicates that the lines l
a
i
,
j and l
b
i
,
j are not parallel and that they therefore intersect each other in one point. The solution w
i
,
j
+ we look for will be at this intersection and the corresponding β
00
p is calculated by replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in the equation ρ
2
(
1
)
(
w
i
,
j
+
)
=
0 using (42). For the case where b
i
,
j
T
n
a
i
,
j
=
0 in (43), which is equivalent to a
i
,
j
∥
b
i
,
j, we conclude that the lines l
a
i
,
j and l
b
i
,
j are parallel. This means that if the solution w
i
,
j
+ lies on the line l
a
i
,
j it necessarily also lies on the line l
b
i
,
j and the parameter β
00
p is calculated by solving problem (28) for the current case and replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) using (42).(ii)
Ifa
i
,
j
=
0
∧
b
i
,
j
≠
0, we have the following.Sincea
i
,
j vanishes, there is no need to consider the line l
a
i
,
j and we focus instead on the line l
b
i
,
j on which the solution w
i
,
j
+ has to lie too. We therefore assume the solution has the form(44)
w
i
,
j
+
=
w
0
i
,
j
+
+
α
2
b
i
,
j
+
β
00
p
n
b
i
,
j
.Replacing this time the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in the equation ρ
2
(
1
)
(
w
i
,
j
+
)
=
0 with the help of (44) results in(45)
α
2
=
-
ρ
2
1
w
0
i
,
j
+
b
i
,
j
2
.After solving problem (28) and replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) using (44) we get (46)
β
00
p
=
n
b
i
,
j
T
-
z
i
,
j
n
b
i
,
j
2
.
(iii)
Ifa
i
,
j
=
0
∧
b
i
,
j
=
0, we have the following.This time botha
i
,
j and b
i
,
j vanish and the corresponding lines l
a
i
,
j and l
b
i
,
j do not play a role in the derivation of the solution w
i
,
j
+ anymore. We can directly solve problem (28) without any prior assumption, similarly as done before in other cases, and we get (47)
w
i
,
j
+
=
w
0
i
,
j
+
-
z
i
,
j
.The derivation of the explicit solutions for cases 18, 19, 20, 24, and 27 in Table 1 is done similarly to case 15 explained above.Case  21. There is now one case left, which was not discussed so far. The corresponding condition from Table 1 is (48)
ρ
1
w
i
,
j
+
=
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
=
0
,
ρ
2
1
w
i
,
j
+
=
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
=
0
,
ρ
2
2
w
i
,
j
+
=
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
=
0
.This time we have three equal signs and we have to distinguish the situations a
i
,
j
≠
0, a
i
,
j
=
0
∧
b
i
,
j
≠
0, a
i
,
j
=
0
∧
b
i
,
j
=
0
∧
c
i
,
j
≠
0, and a
i
,
j
=
0
∧
b
i
,
j
=
0
∧
c
i
,
j
=
0. (i)
Ifa
i
,
j
≠
0, we have the following.Sinceρ
1
(
w
i
,
j
+
)
=
0, ρ
2
(
1
)
(
w
i
,
j
+
)
=
0, and ρ
2
(
2
)
(
w
i
,
j
+
)
=
0, we see that the solution w
i
,
j
+ has to lie this time on all three lines l
a
i
,
j, l
b
i
,
j, and l
c
i
,
j. For a
i
,
j
≠
0 we assume that the explicit solution is of the form(49)
w
i
,
j
+
=
w
0
i
,
j
+
+
α
1
a
i
,
j
+
β
000
n
a
i
,
j
.Like in the paragraphs before, replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in ρ
1
(
w
i
,
j
+
)
=
0 by making use of (49) leads to α
1 being equal to the one in (38). To determine β
000 in (49), we argue as follows: (a)
Ifb
i
,
j
T
n
a
i
,
j
≠
0, which is equivalent to a
i
,
j
∦
b
i
,
j, we use (49) in ρ
2
(
1
)
(
w
i
,
j
+
)
=
0 and get (50)
β
000
=
-
ρ
2
1
w
0
i
,
j
+
-
α
1
b
i
,
j
T
a
i
,
j
b
i
,
j
T
n
a
i
,
j
.
(b)
Ifb
i
,
j
T
n
a
i
,
j
=
0 and c
i
,
j
T
n
a
i
,
j
≠
0, that is, a
i
,
j
∥
b
i
,
j and a
i
,
j
∦
c
i
,
j, we use (49) in ρ
2
(
2
)
(
w
i
,
j
+
)
=
0 and get (51)
β
000
=
-
ρ
2
2
w
0
i
,
j
+
-
α
1
c
i
,
j
T
a
i
,
j
c
i
,
j
T
n
a
i
,
j
.
(c)
Ifb
i
,
j
T
n
a
i
,
j
=
0 and c
i
,
j
T
n
a
i
,
j
=
0, this means that a
i
,
j
∥
b
i
,
j and a
i
,
j
∥
c
i
,
j. Thus all the three lines l
a
i
,
j, l
b
i
,
j, and l
c
i
,
j are parallel and by solving problem (28) and using (49) we get (52)
β
000
=
0
.
(ii)
Ifa
i
,
j
=
0
∧
b
i
,
j
≠
0, we have the following.Similarly as in the discussion of (44) in case 15, we assume the solution is of the form(53)
w
i
,
j
+
=
w
0
i
,
j
+
+
α
2
b
i
,
j
+
β
000
n
b
i
,
j
.Replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in ρ
2
(
1
)
(
w
i
,
j
+
)
=
0 by making use of (53) leads to α
2 being equal to the one in (45). For β
000 in (53), we argue as follows: (a)
Ifc
i
,
j
T
n
b
i
,
j
≠
0, which is equivalent to b
i
,
j
∦
c
i
,
j, we make use of (53) in ρ
2
(
2
)
(
w
i
,
j
+
)
=
0 and get (54)
β
000
=
-
ρ
2
2
w
0
i
,
j
+
-
α
2
c
i
,
j
T
b
i
,
j
c
i
,
j
T
n
b
i
,
j
.
(b)
Ifc
i
,
j
T
n
b
i
,
j
=
0, which means b
i
,
j
∥
c
i
,
j, we get by solving (28) and using (53) (55)
β
000
=
0
.
(iii)
Ifa
i
,
j
=
0
∧
b
i
,
j
=
0
∧
c
i
,
j
≠
0, we have the following.Here,a
i
,
j and b
i
,
j vanish and we focus therefore on the line l
c
i
,
j on which the solution w
i
,
j
+ has to lie. We assume the explicit solution is of the form(56)
w
i
,
j
+
=
w
0
i
,
j
+
+
α
3
c
i
,
j
,where α
3 is calculated by using (56) to replace (
w
i
,
j
+
-
w
0
i
,
j
+
) in ρ
2
(
2
)
(
w
i
,
j
+
)
=
0 and is given by (57)
α
3
=
-
ρ
2
2
w
0
i
,
j
+
c
i
,
j
2
.
(iv)
Ifa
i
,
j
=
0
∧
b
i
,
j
=
0
∧
c
i
,
j
=
0, we have the following.In this situation we do not have to consider the linesl
a
i
,
j, l
b
i
,
j, and l
c
i
,
j at all and solving (28) results in (58)
w
i
,
j
+
=
w
0
i
,
j
+
.


### 3.4.3. Reformulation of the Conditions for Problem (28) with u
~
i
,
j
=
1

Sincew
i
,
j
+ is unknown but w
0
i
,
j
+ in (9) is known, we have to rewrite the conditions for the 27 cases, which are defined by (29), in terms of ρ
1
(
w
0
i
,
j
+
), ρ
2
(
1
)
(
w
0
i
,
j
+
), and ρ
2
(
2
)
(
w
0
i
,
j
+
). This can be done by incorporating the derived explicit solutions and finding proper partitioning of the space spanned by ρ
1
(
w
0
i
,
j
+
), ρ
2
(
1
)
(
w
0
i
,
j
+
), and ρ
2
(
2
)
(
w
0
i
,
j
+
) (see Figure 3). A similar approach was used in the work of Zach et al. [25] to reformulate the conditions for their proposed thresholding scheme.Figure 3
The space spanned byρ
1
w
0
i
,
j
+, ρ
2
1
w
0
i
,
j
+, and ρ
2
2
w
0
i
,
j
+.As one can see in Figure3, we indicated a cube in the coordinate system that will facilitate for us the problem of finding proper partitioning of the considered space. The edge length of the cube is chosen in such a way that the critical values of the reformulated conditions are covered. The reader will understand the meaning of this sentence afterwards in a better way. For later need, we number the faces of the cube as shown in Figure 4.Figure 4
The numbering of the faces for the cube indicated in Figure3.The orientation of the face numbers in Figure4 indicates from which direction we will look at the faces.In the following we want to show how we can rewrite the conditions for the different cases. The same four groups of cases which were discussed in Section3.4.2 are considered here too and for each of them the reformulation of the conditions can be achieved in a similar fashion.Case  1 and Similar Ones. We start again with the first case in Table 1, which belongs to the group of cases that can be handled in the easiest way compared to the other ones. The condition for this case is given in (33) and by using (35) to replace the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in the condition we get (59)
ρ
1
w
i
,
j
+
=
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
,
ρ
2
1
w
i
,
j
+
=
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
,
ρ
2
2
w
i
,
j
+
=
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
,or rewritten(60)
ρ
1
w
0
i
,
j
+
>
-
a
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
ρ
2
1
w
0
i
,
j
+
>
-
b
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
ρ
2
2
w
0
i
,
j
+
>
-
c
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
.One of the critical values we mentioned before concerning the edge length of the cube shown in Figure 3 can be determined now for case number 1, namely, by the values on the right-hand side in (60).The reformulation of the conditions for cases2, 4, 5, 7, 8, 10, and 11 in Table 1 is achieved similarly to this case and the remaining seven critical values, which are necessary to define the edge length of the cube, are determined in the same way as above.Remark 4.
Before jumping to the next group of cases, we want to explain why the cube in Figure3 will help us in finding proper partitioning of the space spanned by ρ
1
(
w
0
i
,
j
+
), ρ
2
(
1
)
(
w
0
i
,
j
+
), and ρ
2
(
2
)
(
w
0
i
,
j
+
). Our aim is to divide the space into 27 parts, each one defining a region for which exactly one of the 27 cases may apply with their already derived explicit solutions. This partitioning should not be done just arbitrarily and has to be meaningful and to conform to the already set up formulations. Before, we mentioned once that we want the eight critical values from the paragraph before to be covered by the cube. The reason for this is that like this it is possible to reduce the partitioning problem of the whole space to a partitioning problem of the cube. The uniquely defined regions of the cube are the ones from the corresponding eight cases discussed in the paragraph above and they define the corners of it. In Figure 5 we show the location of the regions for different case numbers. For simplicity we drew all the regions with the same size, although they can and most likely will have different sizes. Note that the region for case number 21 is not visible in the figure, because it fully lies in the interior of the cube.
Let us focus now on one face of the cube. We know the values for the boundaries of the corner regions from the critical values we got from the previous paragraph. They are fixed and we introduce for them the namesA
v, A
h, B
v, B
h, C
v, C
h, D
v, and D
h as shown in Figure 6.
For face number6, for example, the values D
v and D
h will be the boundary values to region 1 (see Figure 5) and after having a look at (60) we can conclude that D
v
=
-
a
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j and D
h
=
-
b
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j.
For each face it always holds that(61)
A
v
≤
B
v
,
C
v
≤
D
v
,
B
h
≤
D
h
,
A
h
≤
C
h
.For example, to see that C
v
≤
D
v we consider the reformulation of the condition for case number 2:(62)
ρ
1
w
0
i
,
j
+
<
-
a
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
ρ
2
1
w
0
i
,
j
+
>
-
b
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
ρ
2
2
w
0
i
,
j
+
>
-
c
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,which can be derived, as mentioned before, similarly to case number 1. For face number 6 we see in the same manner as before by having a look at Figure 5 and (62) that C
v
=
-
a
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j. Since a
i
,
j
T
x
i
,
j
=
τ
γ
1
a
i
,
j
2
≥
0 (see (31)), we can conclude that C
v
≤
D
v.
Let us consider again an arbitrary face of the cube. Although we have boundaries for the corner regions, the ones for the regions in between are not fully defined. Since a rectangular division of the face seems to be adequate, we investigate the problem on how to achieve such one. In Figure7 a face is depicted with its corner regions and two ways of possible divisions into 9 rectangles are shown, which we call the “clockwise” and “anticlockwise” division method.
Depending on the values forA
v, A
h, B
v, B
h, C
v, C
h, D
v, and D
h sometimes both the clockwise and the anticlockwise division method work, as, for example, in Figure 7. Depending on the situation however, there can be also values for which only the clockwise or the anticlockwise or even none of both methods will work, where in the latter situation another adequate division method has to be defined.
In Table2 we list the different possible situations for the values A
v, A
h, B
v, B
h, C
v, C
h, D
v, and D
h.
Furthermore, we split the last situation S9 into its parts for later use and introduce the corresponding labelling in Table3.
In the following, we go through the different situations given in Tables2 and 3 and decide which kind of method can be used to achieve a division of the face into 9 rectangular regions: (i)
Situation S1 in Table2 is illustrated in Figure 8 for face number 6 and the corresponding region numbers. For this situation, neither the clockwise nor the anticlockwise division method will work, since the corner regions 1 and 5 overlap. An adapted division method, which is depicted in Figure 8, will be used instead for this situation. Note that in this situation certain corner regions, which were defined through the reformulation of the conditions for the group of cases similar to case 1, have to be adapted by removing from them the region where they overlap.
(ii)
For the situations S2, S3, S5, and S6 in Table2 and S9.1, S9.2, S9.3, S9.4, and S9.5 in Table 3 we can use the clockwise division method illustrated in Figure 7. We define the set of situations for which the clockwise division method can be applied by (63)
S
c
≔
S
2
,
S
3
,
S
5
,
S
6
,
S
9.1
,
S
9.2
,
S
9.3
,
S
9.4
,
S
9.5
.
(iii)
Situation S9.7 in Table3 needs again a special treatment to achieve a meaningful division, since neither the clockwise nor the anticlockwise division method can be applied because of the overlapping of the corner regions 2 and 4, if we consider face number 6, for example. The adapted division method, which will be used instead, is shown in Figure 9 for face number 6 and the corresponding region numbers. Note again that also in this situation certain corner regions have to be adapted by removing from them the region where they overlap.
(iv)
Finally, for the situations S4, S7, and S8 in Table2 and S9.6, S9.8, and S9.9 in Table 3 we can use this time the anticlockwise division method illustrated in Figure 7. The set of situations for which the anticlockwise division method can be applied is then given by(64)
S
a
c
≔
S
4
,
S
7
,
S
8
,
S
9.6
,
S
9.8
,
S
9.9
.We are now able to give an idea on how to reformulate the conditions for the remaining group of cases.Table 2
Different possible situations for the valuesA
v, A
h, B
v, B
h, C
v, C
h, D
v, and D
h.
Number
A
v
⋄
D
v
B
v
⋄
C
v
A
h
⋄
D
h
B
h
⋄
C
h
S1
>
>
>
<
S2
>
>
=
<
=
S3
>
>
<
>
<
=
S4
=
>
=
>
<
S5
=
>
=
=
<
=
S6
=
>
=
<
>
<
=
S7
<
>
<
=
>
<
S8
<
>
<
=
=
<
=
S9
<
>
<
=
<
>
<
=Table 3
Splitting of situation S9 in Table2 into its parts.
Number
A
v
⋄
D
v
B
v
⋄
C
v
A
h
⋄
D
h
B
h
⋄
C
h
S9.1
<
>
<
>
S9.2
<
>
<
=
S9.3
<
>
<
<
S9.4
<
=
<
>
S9.5
<
=
<
=
S9.6
<
=
<
<
S9.7
<
<
<
>
S9.8
<
<
<
=
S9.9
<
<
<
<Figure 5
The net of the cube indicated in Figure3 with the corresponding case numbers for the different regions. The face numbers from Figure 4 are visible in the background.Figure 6
One face of the cube with the depicted corner regions and their labelling for the boundary values. The arrows indicate the direction of increasement.Figure 7
A face with its depicted corner regions (a) and the clockwise (b) and anticlockwise (c) division method, respectively.
(a)
(b)
(c)Figure 8
The division method used for situation S1 in Table2. For the illustration we used face number 6 and the corresponding region numbers.Figure 9
The division method used for situation S9.7 in Table3. For the illustration we used again face number 6 and the corresponding region numbers.Case  3 and Similar Ones. The reformulation of the cases in this group with the same approach as for case number 1 will not lead to proper partitioning of the space shown in Figure 3. Instead we will use a more geometrically based approach which depends on the different situations the boundary values can take.First we can derive for sure, independently of which approach we use, that(65)
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
-
y
i
,
j
-
z
i
,
j
≤
a
i
,
j
T
x
i
,
j
.This can be seen by considering Figure 5, from which we can deduce that region 3 has to lie between regions 1 and 2. By using the critical values of these adjacent regions, which can be determined from (60) and (62), we arrive at the inequality above.For the reformulation of the condition with respect toρ
2
(
1
)
(
w
0
i
,
j
+
), we consider the different situations in which the boundary values of the corner regions can be for face number 6.(i)
If the boundary values are in situation S1 given in Table2, the reformulation is given by (66)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
+
x
i
,
j
+
y
i
,
j
-
z
i
,
j
≥
0
.
(ii)
If the boundary values are in one of the situations in the setS
c defined in (63) or in situation S9.7 in Table 3, we get (67)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
.
(iii)
If the boundary values are in one of the situations inS
a
c (64), the reformulation is (68)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
.For the reformulation of the condition with respect to ρ
2
(
2
)
(
w
0
i
,
j
+
), we consider this time the different situations belonging to face number 4.(i)
If the boundary values are in situation S1, the reformulation is given by(69)
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
-
x
i
,
j
-
y
i
,
j
+
z
i
,
j
≥
0
.
(ii)
If the situation of the boundary values is in the setS
c or corresponds to the situation S9.7, we get (70)
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
.
(iii)
Finally, if the situation of the boundary values is inS
a
c, the reformulation is (71)
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
.Finally, the reformulation of the conditions for cases 6, 9, 12, 13, 14, 16, 17, 22, 23, 25, and 26 in Table 1 is achieved similarly to this case.Case  15 and Similar Ones. For this group of cases we will again make use of the different situations of the boundary values, but this time only for face number 6.For the reformulation of the condition with respect toρ
1
(
w
0
i
,
j
+
), we get the following list:(i)
If the boundary values are in situation S1(72)
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
-
z
i
,
j
<
-
a
i
,
j
T
+
x
i
,
j
+
y
i
,
j
.
(ii)
If the situation of the boundary values is in the setS
c
(73)
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
-
z
i
,
j
≤
-
a
i
,
j
T
-
x
i
,
j
+
y
i
,
j
.
(iii)
If the boundary values are in situation S9.7(74)
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
-
z
i
,
j
<
-
a
i
,
j
T
+
x
i
,
j
-
y
i
,
j
.
(iv)
If the situation of the boundary values is in the setS
a
c
(75)
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
-
z
i
,
j
≤
-
a
i
,
j
T
-
x
i
,
j
-
y
i
,
j
.For the reformulation of the condition with respect toρ
2
(
1
)
(
w
0
i
,
j
+
), we get the following list:(i)
If the boundary values are in situation S1(76)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
-
z
i
,
j
<
-
b
i
,
j
T
+
x
i
,
j
+
y
i
,
j
.
(ii)
If the situation of the boundary values is in the setS
c
(77)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
-
z
i
,
j
≤
-
b
i
,
j
T
-
x
i
,
j
-
y
i
,
j
.
(iii)
If the boundary values are in situation S9.7(78)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
-
z
i
,
j
<
-
b
i
,
j
T
-
x
i
,
j
+
y
i
,
j
.
(iv)
If the situation of the boundary values is in the setS
a
c
(79)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
-
z
i
,
j
≤
-
b
i
,
j
T
+
x
i
,
j
-
y
i
,
j
.Finally, we get for the reformulation of the condition with respect to ρ
2
(
2
)
(
w
0
i
,
j
+
) the following list:(i)
If the boundary values are in situation S1(80)
ρ
2
2
w
0
i
,
j
+
+
min
⁡
c
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
c
i
,
j
T
+
x
i
,
j
+
y
i
,
j
-
z
i
,
j
>
0
.
(ii)
If the situation of the boundary values is in the setS
c or S
a
c
(81)
ρ
2
2
w
0
i
,
j
+
+
min
⁡
c
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
c
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
c
i
,
j
T
-
x
i
,
j
+
y
i
,
j
-
z
i
,
j
,
c
i
,
j
T
+
x
i
,
j
+
y
i
,
j
-
z
i
,
j
>
0
.
(iii)
If the boundary values are in situation S9.7(82)
ρ
2
2
w
0
i
,
j
+
+
min
⁡
c
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
c
i
,
j
T
-
x
i
,
j
+
y
i
,
j
-
z
i
,
j
>
0
.The reformulation of the conditions for cases 18, 19, 20, 24, and 27 in Table 1 is done then similarly to case 15 explained above.Case  21. The last case which is left is case number 21. The corresponding region contains the part of the cube which was not covered so far with the corresponding regions of the other cases. Additionally, since it is unfortunately also possible that there are overlappings for some of the regions which appear at different faces, as, for example, the possible overlapping of region 1 with region 11, these overlappings are also assigned to case number 21 to finally achieve proper partitioning.One could get the impression that this handling of the overlappings is mathematically rather grubby, but to legitimate this we recall that an optimal solution of problem (28) will result in rather small values for ρ
1
(
w
i
,
j
+
), ρ
2
(
1
)
(
w
i
,
j
+
), and ρ
2
(
2
)
(
w
i
,
j
+
) and therefore should be ideally close to case 21 in Table 1. Extending therefore the region for case number 21 should not be harmful. Furthermore, for small values of ρ
1
(
w
i
,
j
+
), ρ
2
(
1
)
(
w
i
,
j
+
), and ρ
2
(
2
)
(
w
i
,
j
+
) a reformulation of the condition can lead more likely to overlappings of certain region parts. Therefore it seems to be adequate to assign these overlappings to region 21.It is possible to define region21 in mathematical terms, but since we explained the content of this region already above this will be rather uninspiring and for the implementation an explicit formulation of the region is also not needed, since it can easily be depicted with the help of the other regions.


## 3.4.1. Geometric Interpretation of Problem (28) with u
~
i
,
j
=
1

Before explaining the derivation of the explicit solutions and the needed reformulation of the conditions for the different cases in (29) in more detail, we want to give a deeper insight of the geometric interpretation of problem (28) with u
~
i
,
j
=
1. To facilitate the notation in the following, we introduce the terms(30)
a
i
,
j
=
∇
T
x
+
w
0
+
i
,
j
,
b
i
,
j
=
∂
x
1
x
1
T
x
+
w
0
+
∂
x
2
x
1
T
x
+
w
0
+
i
,
j
,
c
i
,
j
=
∂
x
1
x
2
T
x
+
w
0
+
∂
x
2
x
2
T
x
+
w
0
+
i
,
j
,
(31)
x
i
,
j
=
τ
γ
1
a
i
,
j
,
y
i
,
j
=
τ
γ
2
b
i
,
j
,
z
i
,
j
=
τ
γ
2
c
i
,
j
.Using the definitions in (9) and the notation in (30), we can rewrite ρ
1
(
w
i
,
j
+
), ρ
2
(
1
)
(
w
i
,
j
+
), and ρ
2
(
2
)
(
w
i
,
j
+
) as (32)
ρ
1
w
i
,
j
+
=
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
,
ρ
2
1
w
i
,
j
+
=
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
,
ρ
2
2
w
i
,
j
+
=
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
.Note that with ρ
1
(
w
i
,
j
+
)
=
0 a line l
a
i
,
j is defined with its normal vector being parallel to the vector a
i
,
j. Similarly, ρ
2
(
1
)
(
w
i
,
j
+
)
=
0 defines a line l
b
i
,
j and ρ
2
(
2
)
(
w
i
,
j
+
)
=
0 a line l
c
i
,
j.Now we can argue similar to Zach et al. in [25] for the geometric interpretation. The mathematical structure of their minimization problem for which they derive the thresholding scheme is very similar to the one we have in (28) with u
~
i
,
j
=
1. The first term in (28), u
i
,
j
+
-
u
0
i
,
j
+
2
+
v
i
,
j
+
-
v
0
i
,
j
+
2, is the squared distance of w
i
,
j
+ to w
0
i
,
j
+, and the terms ρ
1
w
i
,
j
+, ρ
2
1
w
i
,
j
+, and ρ
2
2
w
i
,
j
+ define the unsigned distances to the lines l
a
i
,
j, l
b
i
,
j, and l
c
i
,
j, respectively. Considering now all w
δ
i
,
j
+ with a fixed distance δ to w
0
i
,
j
+, we see that problem (28) with u
~
i
,
j
=
1 is minimized for w
δ
i
,
j
+ closest to the three lines l
a
i
,
j, l
b
i
,
j, and l
c
i
,
j. See Figure 2 for an illustration.Figure 2
The geometric interpretation of minimization problem (28) with u
~
i
,
j
=
1.


## 3.4.2. Derivation of the Explicit Solutions for Problem (28) with u
~
i
,
j
=
1

We are ready now to derive for each of the cases an explicit solution by using (28). We show the derivation of the solutions only for four different cases. The explicit solutions for the remaining cases can be calculated in a similar fashion as the presented ones.Case  1 and Similar Ones. Let us consider the first case in Table 1; that is,(33)
ρ
1
w
i
,
j
+
=
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
>
0
,
ρ
2
1
w
i
,
j
+
=
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
>
0
,
ρ
2
2
w
i
,
j
+
=
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
>
0
.Solving problem (28) for this case leads to the equation (34)
1
2
τ
2
w
i
,
j
+
-
w
0
i
,
j
+
+
γ
1
a
i
,
j
+
γ
2
b
i
,
j
+
γ
2
c
i
,
j
=
!
0
,and by using the notation in (31) the explicit solution can then be written as(35)
w
i
,
j
+
=
w
0
i
,
j
+
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
.The derivation of the explicit solutions for cases 2, 4, 5, 7, 8, 10, and 11 in Table 1 is performed similarly. There are however other cases left which need another treatment, as, for example, case number 3.Case  3 and Similar Ones. The condition for case 3 in Table 1 is given by (36)
ρ
1
w
i
,
j
+
=
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
=
0
,
ρ
2
1
w
i
,
j
+
=
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
>
0
,
ρ
2
2
w
i
,
j
+
=
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
>
0
.Here we have to distinguish the two situations a
i
,
j
≠
0 and a
i
,
j
=
0. (i)
Ifa
i
,
j
≠
0, we have the following.Fromρ
1
(
w
i
,
j
+
)
=
0 we see that the solution w
i
,
j
+ has to lie on the line l
a
i
,
j and of course the distance of it to the other lines l
b
i
,
j and l
c
i
,
j should be minimal. Therefore we can assume that the solution is of the form(37)
w
i
,
j
+
=
w
0
i
,
j
+
+
α
1
a
i
,
j
+
β
0
p
p
n
a
i
,
j
,where n
a
i
,
j is the normal to a
i
,
j. Geometrically this means that w
0
i
,
j
+ is first orthogonally projected to the line l
a
i
,
j and afterwards moved along the line such that it minimizes the distance to l
b
i
,
j and l
c
i
,
j. The unknown α
1 is determined by replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in the equation ρ
1
(
w
i
,
j
+
)
=
0 using (37), which leads to(38)
α
1
=
-
ρ
1
w
0
i
,
j
+
a
i
,
j
2
.Solving problem (28) for the current case and replacing again the term (
w
i
,
j
+
-
w
0
i
,
j
+
) using (37) lead after some simple calculations to the second parameter (39)
β
0
p
p
=
n
a
i
,
j
T
-
y
i
,
j
-
z
i
,
j
n
a
i
,
j
2
.
(ii)
Ifa
i
,
j
=
0, we have the following.Sincea
i
,
j vanishes, there is also no line l
a
i
,
j that we have to consider and we can directly derive the solution w
i
,
j
+ by just solving (28), similar to the first case in Table 1, and we get the solution (40)
w
i
,
j
+
=
w
0
i
,
j
+
-
y
i
,
j
-
z
i
,
j
.Finally, we can derive the explicit solution for cases 6, 9, 12, 13, 14, 16, 17, 22, 23, 25, and 26 in Table 1 similar to case 3 above.Case  15 and Similar Ones. Let us now consider case number 15, which was not covered so far. The condition is this time given by (41)
ρ
1
w
i
,
j
+
=
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
=
0
,
ρ
2
1
w
i
,
j
+
=
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
=
0
,
ρ
2
2
w
i
,
j
+
=
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
>
0
,and we see that we have now two equal signs, namely, for ρ
1
(
w
i
,
j
+
) and ρ
2
(
1
)
(
w
i
,
j
+
). This time we have to distinguish the situations a
i
,
j
≠
0, a
i
,
j
=
0
∧
b
i
,
j
≠
0, and a
i
,
j
=
0
∧
b
i
,
j
=
0. (i)
Ifa
i
,
j
≠
0, we have the following.Fromρ
1
(
w
i
,
j
+
)
=
0 and ρ
2
(
1
)
(
w
i
,
j
+
)
=
0 we see that the solution w
i
,
j
+ has to lie on the lines l
a
i
,
j and l
b
i
,
j and should have a minimal distance to the line l
c
i
,
j. Thus, we assume that the explicit solution is of the form(42)
w
i
,
j
+
=
w
0
i
,
j
+
+
α
1
a
i
,
j
+
β
00
p
n
a
i
,
j
.Similar to the last paragraph, replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in the equation ρ
1
(
w
i
,
j
+
)
=
0 by making use of (42) leads to the same α
1 as in (38). For the parameter β
00
p we get(43)
β
00
p
=
-
ρ
2
1
w
0
i
,
j
+
-
α
1
b
i
,
j
T
a
i
,
j
b
i
,
j
T
n
a
i
,
j
,
if
b
i
,
j
T
n
a
i
,
j
≠
0
,
n
a
i
,
j
T
-
z
i
,
j
n
a
i
,
j
2
,
if
b
i
,
j
T
n
a
i
,
j
=
0
.The case where b
i
,
j
T
n
a
i
,
j
≠
0 in (43), which is equivalent to a
i
,
j
∦
b
i
,
j, indicates that the lines l
a
i
,
j and l
b
i
,
j are not parallel and that they therefore intersect each other in one point. The solution w
i
,
j
+ we look for will be at this intersection and the corresponding β
00
p is calculated by replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in the equation ρ
2
(
1
)
(
w
i
,
j
+
)
=
0 using (42). For the case where b
i
,
j
T
n
a
i
,
j
=
0 in (43), which is equivalent to a
i
,
j
∥
b
i
,
j, we conclude that the lines l
a
i
,
j and l
b
i
,
j are parallel. This means that if the solution w
i
,
j
+ lies on the line l
a
i
,
j it necessarily also lies on the line l
b
i
,
j and the parameter β
00
p is calculated by solving problem (28) for the current case and replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) using (42).(ii)
Ifa
i
,
j
=
0
∧
b
i
,
j
≠
0, we have the following.Sincea
i
,
j vanishes, there is no need to consider the line l
a
i
,
j and we focus instead on the line l
b
i
,
j on which the solution w
i
,
j
+ has to lie too. We therefore assume the solution has the form(44)
w
i
,
j
+
=
w
0
i
,
j
+
+
α
2
b
i
,
j
+
β
00
p
n
b
i
,
j
.Replacing this time the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in the equation ρ
2
(
1
)
(
w
i
,
j
+
)
=
0 with the help of (44) results in(45)
α
2
=
-
ρ
2
1
w
0
i
,
j
+
b
i
,
j
2
.After solving problem (28) and replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) using (44) we get (46)
β
00
p
=
n
b
i
,
j
T
-
z
i
,
j
n
b
i
,
j
2
.
(iii)
Ifa
i
,
j
=
0
∧
b
i
,
j
=
0, we have the following.This time botha
i
,
j and b
i
,
j vanish and the corresponding lines l
a
i
,
j and l
b
i
,
j do not play a role in the derivation of the solution w
i
,
j
+ anymore. We can directly solve problem (28) without any prior assumption, similarly as done before in other cases, and we get (47)
w
i
,
j
+
=
w
0
i
,
j
+
-
z
i
,
j
.The derivation of the explicit solutions for cases 18, 19, 20, 24, and 27 in Table 1 is done similarly to case 15 explained above.Case  21. There is now one case left, which was not discussed so far. The corresponding condition from Table 1 is (48)
ρ
1
w
i
,
j
+
=
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
=
0
,
ρ
2
1
w
i
,
j
+
=
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
=
0
,
ρ
2
2
w
i
,
j
+
=
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
w
i
,
j
+
-
w
0
i
,
j
+
=
0
.This time we have three equal signs and we have to distinguish the situations a
i
,
j
≠
0, a
i
,
j
=
0
∧
b
i
,
j
≠
0, a
i
,
j
=
0
∧
b
i
,
j
=
0
∧
c
i
,
j
≠
0, and a
i
,
j
=
0
∧
b
i
,
j
=
0
∧
c
i
,
j
=
0. (i)
Ifa
i
,
j
≠
0, we have the following.Sinceρ
1
(
w
i
,
j
+
)
=
0, ρ
2
(
1
)
(
w
i
,
j
+
)
=
0, and ρ
2
(
2
)
(
w
i
,
j
+
)
=
0, we see that the solution w
i
,
j
+ has to lie this time on all three lines l
a
i
,
j, l
b
i
,
j, and l
c
i
,
j. For a
i
,
j
≠
0 we assume that the explicit solution is of the form(49)
w
i
,
j
+
=
w
0
i
,
j
+
+
α
1
a
i
,
j
+
β
000
n
a
i
,
j
.Like in the paragraphs before, replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in ρ
1
(
w
i
,
j
+
)
=
0 by making use of (49) leads to α
1 being equal to the one in (38). To determine β
000 in (49), we argue as follows: (a)
Ifb
i
,
j
T
n
a
i
,
j
≠
0, which is equivalent to a
i
,
j
∦
b
i
,
j, we use (49) in ρ
2
(
1
)
(
w
i
,
j
+
)
=
0 and get (50)
β
000
=
-
ρ
2
1
w
0
i
,
j
+
-
α
1
b
i
,
j
T
a
i
,
j
b
i
,
j
T
n
a
i
,
j
.
(b)
Ifb
i
,
j
T
n
a
i
,
j
=
0 and c
i
,
j
T
n
a
i
,
j
≠
0, that is, a
i
,
j
∥
b
i
,
j and a
i
,
j
∦
c
i
,
j, we use (49) in ρ
2
(
2
)
(
w
i
,
j
+
)
=
0 and get (51)
β
000
=
-
ρ
2
2
w
0
i
,
j
+
-
α
1
c
i
,
j
T
a
i
,
j
c
i
,
j
T
n
a
i
,
j
.
(c)
Ifb
i
,
j
T
n
a
i
,
j
=
0 and c
i
,
j
T
n
a
i
,
j
=
0, this means that a
i
,
j
∥
b
i
,
j and a
i
,
j
∥
c
i
,
j. Thus all the three lines l
a
i
,
j, l
b
i
,
j, and l
c
i
,
j are parallel and by solving problem (28) and using (49) we get (52)
β
000
=
0
.
(ii)
Ifa
i
,
j
=
0
∧
b
i
,
j
≠
0, we have the following.Similarly as in the discussion of (44) in case 15, we assume the solution is of the form(53)
w
i
,
j
+
=
w
0
i
,
j
+
+
α
2
b
i
,
j
+
β
000
n
b
i
,
j
.Replacing the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in ρ
2
(
1
)
(
w
i
,
j
+
)
=
0 by making use of (53) leads to α
2 being equal to the one in (45). For β
000 in (53), we argue as follows: (a)
Ifc
i
,
j
T
n
b
i
,
j
≠
0, which is equivalent to b
i
,
j
∦
c
i
,
j, we make use of (53) in ρ
2
(
2
)
(
w
i
,
j
+
)
=
0 and get (54)
β
000
=
-
ρ
2
2
w
0
i
,
j
+
-
α
2
c
i
,
j
T
b
i
,
j
c
i
,
j
T
n
b
i
,
j
.
(b)
Ifc
i
,
j
T
n
b
i
,
j
=
0, which means b
i
,
j
∥
c
i
,
j, we get by solving (28) and using (53) (55)
β
000
=
0
.
(iii)
Ifa
i
,
j
=
0
∧
b
i
,
j
=
0
∧
c
i
,
j
≠
0, we have the following.Here,a
i
,
j and b
i
,
j vanish and we focus therefore on the line l
c
i
,
j on which the solution w
i
,
j
+ has to lie. We assume the explicit solution is of the form(56)
w
i
,
j
+
=
w
0
i
,
j
+
+
α
3
c
i
,
j
,where α
3 is calculated by using (56) to replace (
w
i
,
j
+
-
w
0
i
,
j
+
) in ρ
2
(
2
)
(
w
i
,
j
+
)
=
0 and is given by (57)
α
3
=
-
ρ
2
2
w
0
i
,
j
+
c
i
,
j
2
.
(iv)
Ifa
i
,
j
=
0
∧
b
i
,
j
=
0
∧
c
i
,
j
=
0, we have the following.In this situation we do not have to consider the linesl
a
i
,
j, l
b
i
,
j, and l
c
i
,
j at all and solving (28) results in (58)
w
i
,
j
+
=
w
0
i
,
j
+
.


## 3.4.3. Reformulation of the Conditions for Problem (28) with u
~
i
,
j
=
1

Sincew
i
,
j
+ is unknown but w
0
i
,
j
+ in (9) is known, we have to rewrite the conditions for the 27 cases, which are defined by (29), in terms of ρ
1
(
w
0
i
,
j
+
), ρ
2
(
1
)
(
w
0
i
,
j
+
), and ρ
2
(
2
)
(
w
0
i
,
j
+
). This can be done by incorporating the derived explicit solutions and finding proper partitioning of the space spanned by ρ
1
(
w
0
i
,
j
+
), ρ
2
(
1
)
(
w
0
i
,
j
+
), and ρ
2
(
2
)
(
w
0
i
,
j
+
) (see Figure 3). A similar approach was used in the work of Zach et al. [25] to reformulate the conditions for their proposed thresholding scheme.Figure 3
The space spanned byρ
1
w
0
i
,
j
+, ρ
2
1
w
0
i
,
j
+, and ρ
2
2
w
0
i
,
j
+.As one can see in Figure3, we indicated a cube in the coordinate system that will facilitate for us the problem of finding proper partitioning of the considered space. The edge length of the cube is chosen in such a way that the critical values of the reformulated conditions are covered. The reader will understand the meaning of this sentence afterwards in a better way. For later need, we number the faces of the cube as shown in Figure 4.Figure 4
The numbering of the faces for the cube indicated in Figure3.The orientation of the face numbers in Figure4 indicates from which direction we will look at the faces.In the following we want to show how we can rewrite the conditions for the different cases. The same four groups of cases which were discussed in Section3.4.2 are considered here too and for each of them the reformulation of the conditions can be achieved in a similar fashion.Case  1 and Similar Ones. We start again with the first case in Table 1, which belongs to the group of cases that can be handled in the easiest way compared to the other ones. The condition for this case is given in (33) and by using (35) to replace the term (
w
i
,
j
+
-
w
0
i
,
j
+
) in the condition we get (59)
ρ
1
w
i
,
j
+
=
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
,
ρ
2
1
w
i
,
j
+
=
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
,
ρ
2
2
w
i
,
j
+
=
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
,or rewritten(60)
ρ
1
w
0
i
,
j
+
>
-
a
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
ρ
2
1
w
0
i
,
j
+
>
-
b
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
ρ
2
2
w
0
i
,
j
+
>
-
c
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
.One of the critical values we mentioned before concerning the edge length of the cube shown in Figure 3 can be determined now for case number 1, namely, by the values on the right-hand side in (60).The reformulation of the conditions for cases2, 4, 5, 7, 8, 10, and 11 in Table 1 is achieved similarly to this case and the remaining seven critical values, which are necessary to define the edge length of the cube, are determined in the same way as above.Remark 4.
Before jumping to the next group of cases, we want to explain why the cube in Figure3 will help us in finding proper partitioning of the space spanned by ρ
1
(
w
0
i
,
j
+
), ρ
2
(
1
)
(
w
0
i
,
j
+
), and ρ
2
(
2
)
(
w
0
i
,
j
+
). Our aim is to divide the space into 27 parts, each one defining a region for which exactly one of the 27 cases may apply with their already derived explicit solutions. This partitioning should not be done just arbitrarily and has to be meaningful and to conform to the already set up formulations. Before, we mentioned once that we want the eight critical values from the paragraph before to be covered by the cube. The reason for this is that like this it is possible to reduce the partitioning problem of the whole space to a partitioning problem of the cube. The uniquely defined regions of the cube are the ones from the corresponding eight cases discussed in the paragraph above and they define the corners of it. In Figure 5 we show the location of the regions for different case numbers. For simplicity we drew all the regions with the same size, although they can and most likely will have different sizes. Note that the region for case number 21 is not visible in the figure, because it fully lies in the interior of the cube.
Let us focus now on one face of the cube. We know the values for the boundaries of the corner regions from the critical values we got from the previous paragraph. They are fixed and we introduce for them the namesA
v, A
h, B
v, B
h, C
v, C
h, D
v, and D
h as shown in Figure 6.
For face number6, for example, the values D
v and D
h will be the boundary values to region 1 (see Figure 5) and after having a look at (60) we can conclude that D
v
=
-
a
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j and D
h
=
-
b
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j.
For each face it always holds that(61)
A
v
≤
B
v
,
C
v
≤
D
v
,
B
h
≤
D
h
,
A
h
≤
C
h
.For example, to see that C
v
≤
D
v we consider the reformulation of the condition for case number 2:(62)
ρ
1
w
0
i
,
j
+
<
-
a
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
ρ
2
1
w
0
i
,
j
+
>
-
b
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
ρ
2
2
w
0
i
,
j
+
>
-
c
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,which can be derived, as mentioned before, similarly to case number 1. For face number 6 we see in the same manner as before by having a look at Figure 5 and (62) that C
v
=
-
a
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j. Since a
i
,
j
T
x
i
,
j
=
τ
γ
1
a
i
,
j
2
≥
0 (see (31)), we can conclude that C
v
≤
D
v.
Let us consider again an arbitrary face of the cube. Although we have boundaries for the corner regions, the ones for the regions in between are not fully defined. Since a rectangular division of the face seems to be adequate, we investigate the problem on how to achieve such one. In Figure7 a face is depicted with its corner regions and two ways of possible divisions into 9 rectangles are shown, which we call the “clockwise” and “anticlockwise” division method.
Depending on the values forA
v, A
h, B
v, B
h, C
v, C
h, D
v, and D
h sometimes both the clockwise and the anticlockwise division method work, as, for example, in Figure 7. Depending on the situation however, there can be also values for which only the clockwise or the anticlockwise or even none of both methods will work, where in the latter situation another adequate division method has to be defined.
In Table2 we list the different possible situations for the values A
v, A
h, B
v, B
h, C
v, C
h, D
v, and D
h.
Furthermore, we split the last situation S9 into its parts for later use and introduce the corresponding labelling in Table3.
In the following, we go through the different situations given in Tables2 and 3 and decide which kind of method can be used to achieve a division of the face into 9 rectangular regions: (i)
Situation S1 in Table2 is illustrated in Figure 8 for face number 6 and the corresponding region numbers. For this situation, neither the clockwise nor the anticlockwise division method will work, since the corner regions 1 and 5 overlap. An adapted division method, which is depicted in Figure 8, will be used instead for this situation. Note that in this situation certain corner regions, which were defined through the reformulation of the conditions for the group of cases similar to case 1, have to be adapted by removing from them the region where they overlap.
(ii)
For the situations S2, S3, S5, and S6 in Table2 and S9.1, S9.2, S9.3, S9.4, and S9.5 in Table 3 we can use the clockwise division method illustrated in Figure 7. We define the set of situations for which the clockwise division method can be applied by (63)
S
c
≔
S
2
,
S
3
,
S
5
,
S
6
,
S
9.1
,
S
9.2
,
S
9.3
,
S
9.4
,
S
9.5
.
(iii)
Situation S9.7 in Table3 needs again a special treatment to achieve a meaningful division, since neither the clockwise nor the anticlockwise division method can be applied because of the overlapping of the corner regions 2 and 4, if we consider face number 6, for example. The adapted division method, which will be used instead, is shown in Figure 9 for face number 6 and the corresponding region numbers. Note again that also in this situation certain corner regions have to be adapted by removing from them the region where they overlap.
(iv)
Finally, for the situations S4, S7, and S8 in Table2 and S9.6, S9.8, and S9.9 in Table 3 we can use this time the anticlockwise division method illustrated in Figure 7. The set of situations for which the anticlockwise division method can be applied is then given by(64)
S
a
c
≔
S
4
,
S
7
,
S
8
,
S
9.6
,
S
9.8
,
S
9.9
.We are now able to give an idea on how to reformulate the conditions for the remaining group of cases.Table 2
Different possible situations for the valuesA
v, A
h, B
v, B
h, C
v, C
h, D
v, and D
h.
Number
A
v
⋄
D
v
B
v
⋄
C
v
A
h
⋄
D
h
B
h
⋄
C
h
S1
>
>
>
<
S2
>
>
=
<
=
S3
>
>
<
>
<
=
S4
=
>
=
>
<
S5
=
>
=
=
<
=
S6
=
>
=
<
>
<
=
S7
<
>
<
=
>
<
S8
<
>
<
=
=
<
=
S9
<
>
<
=
<
>
<
=Table 3
Splitting of situation S9 in Table2 into its parts.
Number
A
v
⋄
D
v
B
v
⋄
C
v
A
h
⋄
D
h
B
h
⋄
C
h
S9.1
<
>
<
>
S9.2
<
>
<
=
S9.3
<
>
<
<
S9.4
<
=
<
>
S9.5
<
=
<
=
S9.6
<
=
<
<
S9.7
<
<
<
>
S9.8
<
<
<
=
S9.9
<
<
<
<Figure 5
The net of the cube indicated in Figure3 with the corresponding case numbers for the different regions. The face numbers from Figure 4 are visible in the background.Figure 6
One face of the cube with the depicted corner regions and their labelling for the boundary values. The arrows indicate the direction of increasement.Figure 7
A face with its depicted corner regions (a) and the clockwise (b) and anticlockwise (c) division method, respectively.
(a)
(b)
(c)Figure 8
The division method used for situation S1 in Table2. For the illustration we used face number 6 and the corresponding region numbers.Figure 9
The division method used for situation S9.7 in Table3. For the illustration we used again face number 6 and the corresponding region numbers.Case  3 and Similar Ones. The reformulation of the cases in this group with the same approach as for case number 1 will not lead to proper partitioning of the space shown in Figure 3. Instead we will use a more geometrically based approach which depends on the different situations the boundary values can take.First we can derive for sure, independently of which approach we use, that(65)
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
-
y
i
,
j
-
z
i
,
j
≤
a
i
,
j
T
x
i
,
j
.This can be seen by considering Figure 5, from which we can deduce that region 3 has to lie between regions 1 and 2. By using the critical values of these adjacent regions, which can be determined from (60) and (62), we arrive at the inequality above.For the reformulation of the condition with respect toρ
2
(
1
)
(
w
0
i
,
j
+
), we consider the different situations in which the boundary values of the corner regions can be for face number 6.(i)
If the boundary values are in situation S1 given in Table2, the reformulation is given by (66)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
+
x
i
,
j
+
y
i
,
j
-
z
i
,
j
≥
0
.
(ii)
If the boundary values are in one of the situations in the setS
c defined in (63) or in situation S9.7 in Table 3, we get (67)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
.
(iii)
If the boundary values are in one of the situations inS
a
c (64), the reformulation is (68)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
.For the reformulation of the condition with respect to ρ
2
(
2
)
(
w
0
i
,
j
+
), we consider this time the different situations belonging to face number 4.(i)
If the boundary values are in situation S1, the reformulation is given by(69)
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
-
x
i
,
j
-
y
i
,
j
+
z
i
,
j
≥
0
.
(ii)
If the situation of the boundary values is in the setS
c or corresponds to the situation S9.7, we get (70)
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
.
(iii)
Finally, if the situation of the boundary values is inS
a
c, the reformulation is (71)
ρ
2
2
w
0
i
,
j
+
+
c
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
>
0
.Finally, the reformulation of the conditions for cases 6, 9, 12, 13, 14, 16, 17, 22, 23, 25, and 26 in Table 1 is achieved similarly to this case.Case  15 and Similar Ones. For this group of cases we will again make use of the different situations of the boundary values, but this time only for face number 6.For the reformulation of the condition with respect toρ
1
(
w
0
i
,
j
+
), we get the following list:(i)
If the boundary values are in situation S1(72)
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
-
z
i
,
j
<
-
a
i
,
j
T
+
x
i
,
j
+
y
i
,
j
.
(ii)
If the situation of the boundary values is in the setS
c
(73)
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
-
z
i
,
j
≤
-
a
i
,
j
T
-
x
i
,
j
+
y
i
,
j
.
(iii)
If the boundary values are in situation S9.7(74)
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
-
z
i
,
j
<
-
a
i
,
j
T
+
x
i
,
j
-
y
i
,
j
.
(iv)
If the situation of the boundary values is in the setS
a
c
(75)
ρ
1
w
0
i
,
j
+
+
a
i
,
j
T
-
z
i
,
j
≤
-
a
i
,
j
T
-
x
i
,
j
-
y
i
,
j
.For the reformulation of the condition with respect toρ
2
(
1
)
(
w
0
i
,
j
+
), we get the following list:(i)
If the boundary values are in situation S1(76)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
-
z
i
,
j
<
-
b
i
,
j
T
+
x
i
,
j
+
y
i
,
j
.
(ii)
If the situation of the boundary values is in the setS
c
(77)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
-
z
i
,
j
≤
-
b
i
,
j
T
-
x
i
,
j
-
y
i
,
j
.
(iii)
If the boundary values are in situation S9.7(78)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
-
z
i
,
j
<
-
b
i
,
j
T
-
x
i
,
j
+
y
i
,
j
.
(iv)
If the situation of the boundary values is in the setS
a
c
(79)
ρ
2
1
w
0
i
,
j
+
+
b
i
,
j
T
-
z
i
,
j
≤
-
b
i
,
j
T
+
x
i
,
j
-
y
i
,
j
.Finally, we get for the reformulation of the condition with respect to ρ
2
(
2
)
(
w
0
i
,
j
+
) the following list:(i)
If the boundary values are in situation S1(80)
ρ
2
2
w
0
i
,
j
+
+
min
⁡
c
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
c
i
,
j
T
+
x
i
,
j
+
y
i
,
j
-
z
i
,
j
>
0
.
(ii)
If the situation of the boundary values is in the setS
c or S
a
c
(81)
ρ
2
2
w
0
i
,
j
+
+
min
⁡
c
i
,
j
T
-
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
c
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
c
i
,
j
T
-
x
i
,
j
+
y
i
,
j
-
z
i
,
j
,
c
i
,
j
T
+
x
i
,
j
+
y
i
,
j
-
z
i
,
j
>
0
.
(iii)
If the boundary values are in situation S9.7(82)
ρ
2
2
w
0
i
,
j
+
+
min
⁡
c
i
,
j
T
+
x
i
,
j
-
y
i
,
j
-
z
i
,
j
,
c
i
,
j
T
-
x
i
,
j
+
y
i
,
j
-
z
i
,
j
>
0
.The reformulation of the conditions for cases 18, 19, 20, 24, and 27 in Table 1 is done then similarly to case 15 explained above.Case  21. The last case which is left is case number 21. The corresponding region contains the part of the cube which was not covered so far with the corresponding regions of the other cases. Additionally, since it is unfortunately also possible that there are overlappings for some of the regions which appear at different faces, as, for example, the possible overlapping of region 1 with region 11, these overlappings are also assigned to case number 21 to finally achieve proper partitioning.One could get the impression that this handling of the overlappings is mathematically rather grubby, but to legitimate this we recall that an optimal solution of problem (28) will result in rather small values for ρ
1
(
w
i
,
j
+
), ρ
2
(
1
)
(
w
i
,
j
+
), and ρ
2
(
2
)
(
w
i
,
j
+
) and therefore should be ideally close to case 21 in Table 1. Extending therefore the region for case number 21 should not be harmful. Furthermore, for small values of ρ
1
(
w
i
,
j
+
), ρ
2
(
1
)
(
w
i
,
j
+
), and ρ
2
(
2
)
(
w
i
,
j
+
) a reformulation of the condition can lead more likely to overlappings of certain region parts. Therefore it seems to be adequate to assign these overlappings to region 21.It is possible to define region21 in mathematical terms, but since we explained the content of this region already above this will be rather uninspiring and for the implementation an explicit formulation of the region is also not needed, since it can easily be depicted with the help of the other regions.


## 3.5. Resolvent Operators for Problem (12)

This section is very similar to Section3.4 and we mainly have to just replace w
+ by w
- and u
~ by 1
-
u
~. To achieve a smooth extension of w
- to the domain Σ, we consider now instead of problem (12) the following problem:(83)
min
w
-
⁡
∫
Ω
f
w
-
1
-
u
~
x
+
μ
s
w
-
d
x
.The spaces X and Y and their inner products are defined in the same way as in Section 3.4 and the discretization of (83) is performed in the same manner as before.The resolvent operator with respect toF
∗ is identical to (27). Finally, compared to the section before, there are only very slight changes of the resolvent operator with respect to G.


## 4. Implementation and Results

Although we have a convex minimization problem with respect to the motion segmentation functionu
~ in (10) and linearised the fidelity term f in (8), we should keep in mind that the overall minimization problem remains nonconvex and that we have to update w
0
+ and w
0
- regularly. Therefore a coarse-to-fine strategy is applied to avoid the risk of getting stuck in a local minimum during the optimisation. The minimal size of the images at the coarsest level is set to n
m
i
n
=
32 and the scaling factor for the pyramid to ξ
=
0.9. At the coarsest level we initialise the displacement fields w
+ and w
- trivially with 0.Because of Remark1 in Section 2.2 one would expect that the choice of the initialisation for the function u
~ can be arbitrary. Indeed, whatever initialisation we choose, a global minimizer for the motion segmentation can be found forfixed
w
+ and w
-. But in our case the values for w
+ and w
- are updated regularly during optimisation.The final iteration scheme consists of two loops. The outer loop iterates over the pyramid levels. In each level an inner loop updates the values foru
~, w
+, and w
-, following steps (1)–(3) in Section 3.1. Thus, in each iteration of this inner loop one step of the primal-dual Algorithm 1 from [17] is performed with the appropriate resolvent operators from the previous sections to update u
~, w
+, and w
-. To achieve a better convergence of u
~ we update w
+ and w
- only every 50th iteration. The motion segmentation Σ is obtained by choosing η
=
0.5 and setting Σ
=
Σ
(
η
)
≔
x
∈
Ω
∣
u
~
x
≥
η (see Remark 1 in Section 2.2).As soon as the finest level is reached, the inner loop is executed until a certain tolerance or the maximum number of iterations is reached. The final displacement field is then obtained by setting(84)
w
x
=
w
+
x
if
x
∈
Σ
,
w
-
x
if
x
∈
Ω
∖
Σ
.Bicubic interpolation is used to calculate the images T
(
x
+
w
±
) during the iterations and to obtain the final registered image T
(
x
+
w
).In Figure10, we show the registration result of a synthetic example, where a textured circle is moving down diagonally. Although our energy functional is convex with respect to the motion segmentation function u
~, the overall minimization task remains a nonconvex problem. To decrease the influence of the nonconvexity during the minimization procedure, certain terms are linearised and a coarse-to-fine strategy is applied. Experiments for the synthetic example in Figure 10 show that due to this workaround different initialisation for u
~ could be used to achieve similar results, as shown in Figure 11. In order to show the effect of the smoothness parameter ν for the TV norm, we show in Figure 12 that we can obtain smoother results as ν increases.Figure 10
The template (a) and the reference (b) images superimposed with the motion segmentation (green curves) and the displacement fieldw (c). The parameter values used for this example: γ
1
=
4, γ
2
=
4, μ
=
0.5, and ν
=
0.3.
(a)
(b)
(c)Figure 11
The motion segmentation resultsu
~ (b) and the corresponding displacement field results w (c) obtained when using different initialisations for u
~ (a) for the example and parameters in Figure 10. For the results in Figure 10 itself the reference image was used for initialisation.
(a)
(b)
(c)Figure 12
From top to bottom: the effect of the smoothness parameterν when ν
= 0.01, 0.05, and 0.2. Here we show the reference image superimposed with the contour of thresholded motion segmentation u
~ (a), u
~ (b), and the corresponding displacement field w (c) for the example in Figure 10. The other parameters are again γ
1
=
4, γ
2
=
4, and μ
=
0.5.
(a)
(b)
(c)In order to show the performance of the proposed method, we show qualitative and quantitative results. Since we are interested in the discontinuities of the displacement field, we compare the proposed method to methods that are able to preserve discontinuities in the displacement field, in this case, the demon algorithm with anisotropic diffusion filtering [26], the registration algorithm of Brox et al. [8], and our previous work [16].We show the qualitative comparison in Figure13. As clearly seen, the proposed method can achieve more crisp discontinuities in the displacement field compared to the demon algorithm with anisotropic diffusion filtering and the registration algorithm of Brox et al. Furthermore, the proposed method managed nicely to separate the motion of the abdominal wall and the one of the organs.Figure 13
The template image (a), the reference image (b), the difference image (c), and the motion segmentation result (d). The displacement field for the demon algorithm with anisotropic diffusion filtering is shown in (e), the one for the registration algorithm of Brox et al. in (f), the one for our previous method in (g), and finally the one of our method in (h).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)In Figure14 a quantitative evaluation is shown for 22 different liver image pairs from the sequences of [3]. We chose the parameters for all the methods by optimising them with respect to these image pairs. The parameters of our method were set to γ
1
=
4, γ
2
=
1, μ
=
0.2, and ν
=
0.1. For the demon algorithm with anisotropic diffusion filtering we could use the suggested parameters, for Brox et al.’s method we used γ
=
5, α
=
80, and σ
=
0.9, and finally for our previous method we used the parameters suggested there, namely, γ
=
0.4, μ
=
0.05, λ
=
0.05, θ
=
4, and ϵ
=
0.00001.Figure 14
Quantitative evaluation for22 pairs of liver images with a discontinuous displacement field. Comparison of the MSE (a) and NMI values (b).
(a)
(b)To quantitatively assess the performance of the four methods, we calculated the mean squared error (MSE) and the normalised mutual information (NMI), with the grey values scaled from0 to 1. For all our examples the proposed method performed better than the demon algorithm with anisotropic diffusion and the registration algorithm of Brox et al. We used the Kolmogorov-Smirnov test to check for normality of the results using the R Software package (Version 2.10.1). We considered a significance level of 5% as significant. The t-test showed that the proposed method delivered significantly better results than the demon algorithm with anisotropic diffusion filtering and the method of Brox et al. with both p
<
0.05. There was however no significant difference in the performance of the proposed method compared to our previous method [16] for the MSE and NMI.In a next step we compared the running times of the proposed method to the ones of our previous work [16] for the 22 liver image pairs. Both methods were implemented in MATLAB and the experiments were performed on a 64-bit Linux system with 1.2 GHz. We used again the parameters γ
1
=
4, γ
2
=
1, μ
=
0.2, and ν
=
0.1 for the proposed method and for our previous work the parameters γ
=
0.4, μ
=
0.05, λ
=
0.05, θ
=
4, and ϵ
=
0.00001. In [16] the displacement fields w
+ and w
- were updated every 10th iteration. To provide a fair comparison, we therefore used the same update frequency for the proposed new method. Furthermore, the maximum number of iterations in the finest level is set to 30000 for both methods. The timing results for the 22 liver image pairs are shown in Table 4 together with the mean and standard deviation. The running times of both methods are comparable and the proposed new method performed with around 100 s slightly faster than the old method. The old method makes use of some optimised in-built MATLAB routines, whereas for the proposed method there is still a high potential to optimise the code.Table 4
Comparison of the running times of the proposed method to the ones of our previous work [16] for 22 liver image pairs. The times are given in seconds.
Ex. number
Previous method [16]
Proposed new method
1
1055.13
959.33
2
1053.55
951.41
3
1060.33
948.92
4
1059.57
996.78
5
1052.58
951.85
6
1064.84
955.43
7
1070.07
990.57
8
1060.53
944.11
9
1082.90
949.47
10
1079.99
949.23
11
1073.25
954.15
12
1048.56
961.26
13
1059.35
956.23
14
1061.09
949.74
15
1065.24
952.35
16
1072.26
949.80
17
1075.99
980.85
18
1068.88
951.10
19
1071.08
985.62
20
1064.42
949.60
21
1057.28
946.05
22
1060.19
979.33
Mean
1064.41
959.69
Std.
9.03
15.76


## 5. Conclusion

In this paper we presented a primal-dual method for discontinuity preserving nonrigid registration that makes use of the continuous cuts framework. The so-gained motion segmentation influences the motion estimation positively by sharpening the discontinuities in the displacement field. The minimization of the energy functional was implemented in a coarse-to-fine strategy and exploits the rapidity of the primal-dual algorithm studied in [17]. The experimental results demonstrated desirable performance of the proposed method in comparison with those of the demon algorithm with anisotropic diffusion filtering [26] and the registration algorithm of Brox et al. [8].Large displacements of deforming organs can cause misregistrations even when using a coarse-to-fine approach. Brox and Malik proposed in a very recent publication [27] to include point correspondences from descriptor matching into the variational optical flow formulation. In future work we plan on similarly including known corresponding landmarks into the functional as hard constraints as we have preliminarily done in [28]. However, we plan on using ourTracking the Invisible approach [29] for locating the matching point correspondences.



---

*Source: 9504949-2016-09-19.xml*