Noname manuscript No.

(will be inserted by the editor)

Variational Methods for Normal Integration

Yvain Quéau · Jean-Denis Durou · Jean-François Aujol

the date of receipt and acceptance should be inserted later

Abstract The need for an efficient method of integra- to estimate the values of a function z : R2 → R, over a
tion of a dense normal field is inspired by several com- set Ω ⊂ R2 (reconstruction domain) where an estimate
puter vision tasks, such as shape-from-shading, pho- g = [p, q]> : Ω → R2 of its gradient ∇z is available.
tometric stereo, deflectometry, etc. Inspired by edge- Formally, we want to solve the following equation in the
preserving methods from image processing, we study in unknown depth map z:
this paper several variational approaches for normal in- >
tegration, with a focus on non-rectangular domains, free ∇z(u, v) = [p(u, v), q(u, v)] , ∀(u, v) ∈ Ω (1)
| {z }
boundary and depth discontinuities. We first introduce g(u,v)
a new discretization for quadratic integration, which is
designed to ensure both fast recovery and the ability In a companion survey paper [48], we have shown
to handle non-rectangular domains with a free bound- that an ideal numerical tool for solving Equation (1)
ary. Yet, with this solver, discontinuous surfaces can be should satisfy the following properties, appart accuracy:
handled only if the scene is first segmented into pieces • PFast : be as fast as possible;
without discontinuity. Hence, we then discuss several • PRobust : be robust to a noisy gradient field;
discontinuity-preserving strategies. Those inspired, re- • PFreeB : be able to handle a free boundary;
spectively, by the Mumford-Shah segmentation method • PDisc : preserve the depth discontinuities;
and by anisotropic diffusion, are shown to be the most • PNoRect : be able to work on a non-rectangular do-
effective for recovering discontinuities. main Ω;
• PNoPar : have no critical parameter to tune.
Keywords 3D-reconstruction, integration, normal
field, gradient field, variational methods, photometric
Contributions. This paper builds upon the previous con-
stereo, shape-from-shading.
ference papers [19, 20, 47] to clarify the building blocks
of variational approaches to the integration problem,
with a view to meeting the largest subset of these re-
1 Introduction quirements. As discussed in Section 2, the variational
framework is well-adapted to this task, thanks to its
In this paper, we study several methods for numerical flexibility. However, these properties are difficult, if not
integration of a gradient field over a 2D-grid. Our aim is impossible, to satisfy simultaneously. In particular, PDisc
Y. Quéau
seems hardly compatible with PFast and PNoPar .
Technical University Munich, Germany Therefore, we first focus in Section 3 on the prop-
E-mail: [Link]@[Link] erties PFreeB and PNoRect . A new discretization strat-
J.-D. Durou egy for normal integration is presented, which is in-
IRIT, Université de Toulouse, France dependent from the shape of the domain and assumes
J.-F. Aujol no particular boundary condition. When used within a
IMB, Université de Bordeaux, Talence, France quadratic variational approach, this discretization strat-
Institut Universitaire de France egy allows to ensure all the desired properties except
2 Yvain Quéau et al.

PDisc . In particular, the numerical solution comes down • F(z) is a fidelity term penalizing the difference be-
to solving a symmetric, diagonally dominant linear sys- tween a corrupted image z 0 and the restored image:
tem, which can be achieved very efficiently using pre-
conditioning techniques. In comparison with our previ-
ZZ
Φ z(u, v) − z 0 (u, v) du dv


F(z) = (3)
ous work [20] which considered only forward finite dif-
ferences and standard Jacobi iterations, the properties (u,v)∈Ω

PRobust and PFast are better satisfied. • R(z) is a regularization term, which usually penal-
In Section 4, we focus more specifically on the inte- izes the gradient of the restored image:
gration problem in the presence of discontinuities. Sev- ZZ
eral variations of well-known models from image pro- R(z) = λ(u, v) Ψ (k∇z(u, v)k) du dv (4)
cessing are empirically compared, while suggesting for (u,v)∈Ω
each of them the appropriate state-of-the-art minimiza-
In (3), Φ is chosen accordingly to the type of corrup-
tion method. Besides the approaches based on total
tion the original image z 0 is affected by. For instance,
variation and non-convex regularization, which we al-
ΦL2 (s) = s2 is the natural choice in the presence of ad-
ready presented, respectively, in [47] and [19], two new
ditive, zero-mean, Gaussian noise, while ΦL1 (s) = |s|
methods inspired by the Mumford-Shah segmentation
can be used in the presence of bi-exponential (Lapla-
method and by anisotropic diffusion are introduced.
cian) noise, which is a rather good model when outliers
They are shown to be particularly effective for handling
come into play (e.g., “salt & pepper” noise).
PDisc , although PFast and PNoPar are lost.
In (4), λ ≥ 0 is a field of weights which control the
These variational methods for normal integration
respective influence of the fidelity and the regulariza-
are based on the same variational framework, which is
tion terms. It can be either manually tuned beforehand
detailed in the next section.
(if λ(u, v) ≡ λ, λ can be seen as a “hyper-parameter”),
or defined as a function of k∇z(u, v)k.
The choice of Ψ must be made accordingly to a de-
2 From Variational Image Restoration to
sired smoothness of the restored image. The quadratic
Variational Normal Integration
penalty ΨL2 (s) = s2 will produce “smooth” images,
while piecewise-constant images are obtained when choos-
In view of the PRobust property, variational methods,
ing the sparsity penalty ΨL0 (s) = 1 − δ(s), with δ(s) =
which aim at estimating the surface by minimization of
1 if s = 0 and δ(s) = 0 otherwise. The latter ap-
a well-chosen criterion, are particularly suited for the
proach preserves the edges, but the numerical solving
integration problem. Hence, we choose the variational
is much more difficult, since the regularization term is
framework as basis for the design of new methods. This
non-smooth and non-convex. Hence, several choices of
choice is also motivated by the fact that the property
regularizers “inbetween” the quadratic and the sparsity
which is the most difficult to ensure is probably PDisc .
ones have been suggested.
Numerous variational methods have been designed for
For instance, the total variation (TV) regularizer
edge-preserving image processing: such methods may
is obtained by setting Ψ(s) = |s|. Efficient numerical
thus be a natural source of inspiration for designing
methods exist for solving this non-smooth, yet convex,
discontinuity-preserving integration methods.
problem. Examples include primal-dual methods [13],
augmented Lagrangian approaches [23], and forward-
backward splittings [40]. The latter can also be adapted
2.1 Variational Methods in Image Processing
to the case where the regularizer Ψ is non-convex, but
For a comprehensive introduction to this literature, we smooth [41]. Such non-convex regularization terms were
refer the reader to [4] and to pioneering papers such shown to be particularly effective for edge-preserving
as [11, 16, 34, 38]. Basically, the idea in edge-preserving image restoration [22, 36, 38].
image restoration is that edges need to be processed in a Another strategy is to stick to quadratic regulariza-
particular way. This is usually achieved by choosing an tion (Ψ = ΨL2 ), but apply it in a non-uniform manner
appropriate energy to minimize, formulating the inverse by tuning the field of weights λ in (4). For instance, set-
problem as the recovery of a restored image z : Ω ⊂ ting λ(u, v) in (4) inversely proportional to k∇z(u, v)k
R2 → R minimizing the energy: yields the “anisotropic diffusion” model by Perona and
Malik [44]. The discontinuity set K can also be auto-
E(z) = F(z) + R(z) (2) matically estimated and discarded by setting λ(u, v) ≡
0 over K and λ(u, v) ≡ λ over Ω\K, in the spirit of
where: Mumford and Shah’s segmentation method [37].
Variational Methods for Normal Integration 3

2.2 Notations We will discuss methods seeking the depth z as the

minimizer of an energy E(z) in the form (2), but with
Although we chose for simplicity to write the variational different choices for F(z) and R(z):
problems in a continuous form, we are overall interested • F(z) now represents a fidelity term penalizing the
in solving discrete problems. Two different discretiza- difference between the gradient of the recovered depth
tion strategies exist. The first one consists in using vari- map z and the datum g:
ational calculus to derive the (continuous) necessary ZZ
optimality condition, then discretize it by finite differ- F(z) = Φ (k∇z(u, v) − g(u, v)k) du dv (5)
ences, and eventually solve the discretized optimality (u,v)∈Ω
condition. The alternative method is to discretize the
• The regularization term R(z) now represents prior
functional itself by finite differences, before solving the
knowledge of the depth1 :
optimality condition associated to the discrete problem. ZZ
2
λ(u, v) z(u, v) − z 0 (u, v)

As shown in [20], the latter approach eases the han- R(z) = (6)
dling of the boundary of Ω, hence we use it as dis- (u,v)∈Ω
cretization strategy. The variational models hereafter
will be presented using the continuous notations, be- where z 0 is the prior, and λ(u, v) ≥ 0 is a user-
cause we find them more readable. The discrete nota- defined, spatially-varying, regularization weight. In
tions will be used only when presenting the numerical this work, we consider for simplicity only the case
solving. Yet, to avoid confusion, we will use caligraphic where λ does not depend on z.
letters for the continuous energies (e.g., E), and capital
letters for their discrete counterparts (e.g., E). With 2.4 Choosing λ and z 0
these conventions, it should be clear whether an opti-
mization problem is discrete or continuous. Hence, we The main purpose of the regularization term R defined
>
will use the same notation ∇z = [∂u z, ∂v z] both for in (6) is to avoid numerical instabilities which may arise
the gradient of z and its finite differences approxima- when considering solely the fidelity term (5): this fi-
tion. delity term depends only on ∇z, and not on z itself,
hence the minimizer of (5) can be estimated only up to
an additive ambiguity.
Besides, one may also want to impose one or sev-
2.3 Proposed Variational Framework eral control points on the surface [30, 33]. This can be
achieved very simply within the proposed variational
In this work, we show how to adapt the aforemen- framework, by setting λ(u, v) ≡ 0 everywhere, except
tioned variational models, originally designed for im- on the control points locations (u, v) where a high value
age restoration, to normal integration. Although both for λ(u, v) must be set and the value z 0 (u, v) is fixed.
these inverse problems are formally very similar, they Another typical situation is when, given both a coarse
are somehow different, for the following reasons: depth estimate and an accurate normal estimate, one
would like to “merge” them in order to create a high-
• The concept of edges in an image to restore is re- quality depth map. Such a problem arises, for instance,
placed by those of depth discontinuities and kinks. when refining the depth map of an RGB-D sensor (e.g.,
• Unlike image processing functionals, our data con- a Kinect) by means of shape-from-shading [42], pho-
sist in an estimate g of the gradient of the unknown tometric stereo [25] or shape-from-polarization [32]. In
z, in lieu of a corrupted version z 0 of z. Therefore, such cases, we may set z 0 to the coarse depth map, and
the fidelity term F(z) will apply to the difference tune λ so as to merge the coarse and fine estimates in
between ∇z and g, and it is the choice of this term the best way. Non-uniform weights may be used, in or-
which will or not allow depth discontinuities. der to lower the influence of outliers in the coarse depth
• Regularization terms are optional here: all the meth- map [25].
ods we discuss basically work even with R(z) ≡ 0, Eventually, in the absence of such priors, we will
but we may use this regularization term to allow in- use the regularization term only to fix the integration
troducing, if available, a prior on the surface (e.g., constant: this is easily achieved by setting an arbitrary
user-defined control points [30, 33] or a rough depth prior (e.g., z 0 (u, v) ≡ 0), along with a small value for λ
estimate obtained using a low-resolution depth sen- (typically, λ(u, v) ≡ λ = 10−6 ).
sor [32]). Such feature “is appreciable, although not 1
We consider only quadratic regularization terms: studying
required” [48]. more robust ones (e.g., L1 norm) is left as perspective.
4 Yvain Quéau et al.

3 Smooth Surfaces For now, let us use a Gaussian approximation for the
noise contained in g6 , i.e., let us assume in the rest of
We first tackle the problem of recovering a “smooth” this section that each datum gu,v , (u, v) ∈ Ω, is equal
depth map z from a noisy estimate g of ∇z. To this to the gradient ∇z(u, v) of the unknown depth map
end, we consider the quadratic variational problem: z, taken at point (u, v), up to a zero-mean additive,
ZZ homoskedastic (same variance at each location (u, v)),
2
min k∇z(u, v) − g(u, v)k Gaussian noise:
z
2
(u,v)∈Ω + λ(u, v) z(u, v) − z 0 (u, v) du dv

(7) gu,v = ∇z(u, v) + (u, v) (10)
When λ ≡ 0, Problem (7) comes down to Horn and   2 
> σ 0
Brook’s model [29]. In that particular case, an infinity where (u, v) ∼ N [0, 0] , and σ is unknown7 .
0 σ2
of solutions z ∈ W 1,2 (Ω) exist, and they differ by an ad-
Now, we need to give a discrete interpretation of the
ditive constant2 . On the other hand, the regularization
gradient operator in (10), through finite differences.
term allows us to guarantee uniqueness of the solution
In order to obtain a second-order accurate discretiza-
as soon as λ is strictly positive almost everywhere3,4 .
tion, we combine forward and backward first-order fi-
If the depth map z is further assumed to be twice
nite differences, i.e. we consider that each measure of
differentiable, the necessary optimality condition as- >
the gradient gu,v = [pu,v , qu,v ] provides us with up
sociated to the continuous optimization problem (7)
to four independent and identically distributed (i.i.d.)
(Euler-Lagrange equation) is written:
statistical observations, depending on the neighborhood
− ∆z + λz = −∇ · g + λz 0 over Ω (8) of (u, v). Indeed, its first component pu,v can be under-
stood either in terms of both forward or backward finite
(∇z − g) · η = 0 over ∂Ω (9)
differences (when both the bottom and the top8 neigh-
with η a normal vector to the boundary ∂Ω of Ω, ∆ bors are inside Ω), by one of both these discretizations
the Laplacian operator, and ∇· the divergence oper- (only one neighbor inside Ω), or by none of these finite
ator. This condition is a linear PDE in z which can differences (no neighbor inside Ω). Formally, we model
be discretized using finite differences. Yet, providing a the p-observations in the following way:
consistent discretization on the boundary of Ω is not
∂ + zu,v
straightforward [26], especially when dealing with non- z u}| {
rectangular domains Ω where many cases have to be pu,v = zu+1,v − zu,v ++
u (u, v),

considered [6]. Hence, we follow a different track, based ∀(u, v) ∈ {(u, v) ∈ Ω | (u + 1, v) ∈ Ω} (11)
| {z }
on the discretization of the functional itself. +
Ωu
−
∂u zu,v

pu,v = zu,v − zu−1,v +−

z }| {
u (u, v),
3.1 Discretizing the Functional
∀(u, v) ∈ {(u, v) ∈ Ω | (u − 1, v) ∈ Ω} (12)
2
| {z }
Instead of a continuous gradient field g : Ω → R over −
Ωu
an open set Ω, we are actually given a finite set of val-
+/−
ues {gu,v = [pu,v , qu,v ]> , (u, v) ∈ Ω}, where the (u, v) where u ∼ N (0, σ 2 ). Hence, rather than considering
represent the pixels of a discrete subset Ω of a regular that we are given |Ω| observations p, our discretization
square 2D-grid5 . Solving the discrete integration prob- handles these data as |Ωu+ |+|Ωu− | observations, some of
lem requires estimating a finite set of values, i.e. the |Ω| them being interpreted in terms of forward differences,
unknown depth values zu,v , (u, v) ∈ Ω (| · | denotes the some in terms of backward differences, some in terms
cardinality), which are stacked columnwise in a vector of both forward and backward differences, the points
z ∈ R|Ω| . without any neighbor in the u-direction being excluded.
2 6
Proof: by developing the terms inside the integral in (7), In 3D-reconstruction applications such as photometric
and integrating by parts, Theorem 6.2.5 in [3] applies with stereo [55], the assumption on the noise should rather be for-
f := −∇ · g and g := g · η. mulated on the images. This will be discussed in more details
3
Proof: by developing the terms inside the integral in (7) in Subsection 4.4.
7
and integrating by parts, Theorem 6.2.2-(ii) in [3] applies with The assumptions of equal variance σ 2 for both compo-
f := −∇ · g + λz 0 and g := g · η. nents and of a diagonal covariance matrix are introduced only
4
This condition makes the matrix of the associated discrete for consistency with the least-squares problem (7). They are
problem strictly diagonally dominant, see Section 3.2. discussed with more care in Subsection 4.4.
5 8
To ease the comparison between the variational and the The u-axis points “downwards”, the v-axis points “to the
discrete problems, we will use the same notation Ω for both right” and the z-axis points from the surface to the camera,
the open set of R2 and the discrete subset of the grid. see Figure 1.
Variational Methods for Normal Integration 5

Symmetrically, the second component q of g corre- instance, the i-th line of D+ u reads:
sponds either to two, one or zero observations:
D+


u i,· =
∂ + zu,v
h i
z v}| {  −1 , |{z}
 0, . . . , 0, |{z} 1 , 0, . . . , 0 if m(i) ∈ Ωu+
qu,v = zu,v+1 − zu,v ++
v (u, v), Position i Position i+1
 >
∀(u, v) ∈ {(u, v) ∈ Ω | (u, v + 1) ∈ Ω} (13)  0 otherwise
| {z }
Ωv+ (17)
∂v− zu,v
where m is the mapping associating linear indices i with
qu,v = zu,v − zu,v−1 +−
z }| {
v (u, v), the pixel coordinates (u, v):
∀(u, v) ∈ {(u, v) ∈ Ω | (u, v − 1) ∈ Ω} (14)
m : {1, . . . , |Ω|} → Ω
(18)
| {z }
Ωv− i 7→ m(i) = (u, v)
+/− Once these matrices are defined, (16) is equal to:
where v ∼ N (0, σ 2 ). Given the Gaussianity of the
+/−
noises u/v , their independence, and the fact that they 1 + 2 2

FL2 (z) = Du z − p + D− uz−p
all share the same standard deviation σ and mean 0, the 2
1 + 2 2

joint likelihood of the observed gradients {gu,v }(u,v) is: + Dv z − q + D− vz−q
2 
L({gu,v , (u, v) ∈ Ω} | {zu,v , (u, v) ∈ Ω}) 1 X X X X
−  pu,v 2 + pu,v 2 
( 2
) 2
1 [∂u+ zu,v − pu,v ] + −
(u,v)∈Ω\Ωu (u,v)∈Ω\Ωu
Y
= √ exp − 
+ 2πσ 2 2σ 2 1  XX XX
(u,v)∈Ωu ( 2
) − qu,v 2 + qu,v 2  (19)
Y 1 [∂u− zu,v − pu,v ] 2
× √ exp − +
(u,v)∈Ω\Ωv −
(u,v)∈Ω\Ωv
− 2πσ 2 2σ 2
(u,v)∈Ωu ( ) The terms in both the last rows of (19) being indepen-
2
Y 1 [∂v+ zu,v − qu,v ] dent from the z-values, they do not influence the actual
× √ exp −
2πσ 2 2σ 2 minimization and will thus be omitted from now on.
(u,v)∈Ωv+ ( 2
) The regularization term (6) is discretized as:
Y 1 [∂v− zu,v − qu,v ]
× √ exp − (15) XX
0
2
2
= Λ z − z0

− 2πσ 2 2σ 2 R(z) = λu,v zu,v − zu,v (20)
(u,v)∈Ωv
(u,v)∈Ω

and hence the maximum-likelihood estimate for the depth with Λ a |Ω|×|Ω| diagonal matrix containing the values
p
values is obtained by minimizing: λu,v , (u, v) ∈ Ω.
Putting it altogether, our quadratic integration method
1 XX  + 2 reads as the minimization of the discrete functional:
FL2 (z) = ∂u zu,v − pu,v
2 +
1 + 2 − 2

(u,v)∈Ωu ! E L2
(z) = D u z − p + D u z − p
XX  2 2
+ ∂u− zu,v − pu,v 1 + 2 2

2
−
−
(u,v)∈Ωu + D v z − q + Dv z − q + Λ z − z0 (21)
2
1 XX  + 2
+ ∂v zu,v − qu,v
2 3.2 Numerical Solution
(u,v)∈Ωv+ !
X X 
−
2
+ ∂v zu,v − qu,v (16)
The optimality condition associated with the discrete
−
(u,v)∈Ωv
functional (21) is a linear equation in z:
where the 21 coefficients are meant to ease the continu- Az = b (22)
ous interpretation: the integral of the fidelity term in (7) where A is a |Ω| × |Ω| symmetric matrix9 :
is approximated by FL2 (z), expressed in (16) as the
L
mean of the forward and the backward discretizations. z }| {
To obtain a more concise representation of this fi- 1 h +> + −> − +> + −> −
i
A = Du Du + Du Du + Dv Dv + Dv Dv
delity term, let us stack the data in two vectors p ∈ R|Ω| 2
and q ∈ R|Ω| . In addition, let us introduce four |Ω|×|Ω| + Λ2 (23)
− −
differentiation matrices D+ +
u , Du , Dv and Dv , asso- 9
A and b are purposely divided by two in order to ease
+/−
ciated with the finite differences operators ∂u/v . For the continuous interpretation of Subsection 3.3.
6 Yvain Quéau et al.

and b is a |Ω| × 1 vector: based on the preconditioned conjugate gradient algo-

rithm. Eventually, let us recall that tuning λ and/or
D D
z }|u { z }|v { manually fixing the values of the prior z0 is necessary
1 h +> − >
i 1 h
+ > −>
i
only to introduce a prior, but not in general. Hence,
b = Du + Du p + Dv + Dv q
2 2 PNoPar is also enforced. In conclusion, all the desired
+ Λ2 z0 (24) properties are satisfied, except PDisc . Let us now pro-
vide additional remarks on the connections between the
The matrix A is sparse: it contains at most five non- proposed discrete approach and a fully variational one.
zero entries per row. In addition, it is diagonal domi-
nant: if (Λ)i,i = 0, the value (A)i,i of a diagonal entry
is equal to the opposite of the sum of the other en- 3.3 Continuous Interpretation
tries (A)i,j , i 6= j, from the same row i. It becomes
strictly superior as soon as (Λ)i,i is strictly positive. System (22) is nothing else than a discrete analogue of
Let us also remark that, when Ω describes a rectangu- the continuous optimality conditions (8) and (9):
lar domain and the regularization weights are uniform Λ2 z = Du p + Dv q + Λ
Lz + |{z} 2 0
z (26)
(λ(u, v) ≡ λ), A is a Toeplitz matrix. Yet, this struc- |{z}
≈−∆z ≈λz
| {z } | {z }
≈−∇·g ≈λz 0
ture is lost in the general case where it can only be
said that A is a sparse, symmetric, diagonal dominant where the matrix-vector products are easily interpreted
(SDD) matrix with at most 5|Ω| non-zero elements. It in terms of the differential operators in the continuous
is positive semi-definite when Λ = 0, and positive defi- formula (8). One major advantage when reasoning from
nite as soon as one of the λu,v is non-zero. the beginning in the discrete setting is that one does not
System (22) can be solved by means of the conju- need to find out how to discretize the natural 12 bound-
gate gradient algorithm. Initialization will not influence ary condition (9), which was already emphasized in [20,
the actual solution, but it may influence the number of 26]. Yet, the identifications in (26) show that both the
iterations required to reach convergence. In our exper- discrete and continuous approaches are equivalent, pro-
iments, we used z 0 as initial guess, yet more elaborate vided that an appropriate discretization of the contin-
initialization strategies may yield faster convergence [6]. uous optimality condition is used. It is thus possible to
To ensure PFast , we used the multigrid preconditioning derive O(5n log(n) log(1/)) algorithms based on the
technique [35], which has a negligible cost of compu- discretization of the Euler-Lagrange equation, contrar-
tation and still bounds the computational complexity ily to what is stated in [26]. The real drawback of such
required to reach a  relative accuracy10 by: approaches does not lie in complexity, but in the dif-
ficult discretization of the boundary condition. This is
O (5n log(n) log(1/)) (25) further explored in the next subsection.

where n = |Ω|11 . This complexity is inbetween the com-

plexities of the approaches based on Sylvester equa- 3.4 Example
tions [26] (O(n1.5 )) and on DCT [53] (O(n log(n))). Be-
sides, these competing methods explicitly require that To clarify the proposed discretization of the integration
Ω is rectangular, while ours does not. problem, let us consider a non-rectangular domain Ω
By construction, the integration method consisting inside a 3 × 3 grid, like the one depicted in Figure 1.
in minimizing (21) satisfies the PRobust property (it The vectorized unknown depth z and the vectorized
is the maximum-likelihood estimate in the presence of components p and q of the gradient write in this case:
zero-mean Gaussian noise). The discretization we in- 
z1,1
 
p1,1
 
q1,1

troduced does not assume any particular shape for Ω, z2,1  p2,1  q2,1 
neither treats the boundary in a specific manner, hence  
z3,1 

p3,1 
  
q3,1 
PFreeB and PNoRect are also satisfied. We also showed  
z1,2 

p1,2 
  
q1,2 
that PFast could be satisfied, using a solving method z= 
 p = p2,2  q = q2,2 
     (27)
z2,2     
10
z3,2  p3,2  q3,2 
In our experiments, the threshold of the stopping criterion      
z1,3  p1,3  q1,3 
is set to  = 10−4 .
11
In (25), the factor 5n is nothing else than the number z2,3 p2,3 q2,3
of non-zero elements in A. Therefore, exploiting sparsity is
12
not as “fruitless” as argued in [26] when it comes to solving As stated in [26], homogeneous Neumann boundary con-
large linear systems faster than using Gaussian elimination ditions of the type ∇z · η = 0, used e.g. in [1], should be
(complexity O(n3 )). avoided.
Variational Methods for Normal Integration 7

v One can observe that this matrix describes the connec-

tivity of the graph representing the discrete domain Ω:
the diagonal elements (L)i,i are the numbers of neigh-
u (1, 1) (1, 2) (1, 3) bors connected to the i-th point, and the off-diagonals
elements (L)i,j are worth −1 if the i-th and j-th points
are connected, 0 otherwise.
(2, 1) (2, 2) (2, 3) Eventually, the matrices Du and Dv defined in (24)
are equal to:

−1 −1 0 0 0 0 0 0
 
(3, 1) (3, 2)
 1 0 −1 0 0 0 0 0 
 
0 1 1 0 0 0 0 0
 
1 0 0 0 −1 −1 0 0 0 
Du =   (34)
2  0 0 0 1 0 −1 0 0  
0 0 0 0 1 1 0 0
 
Fig. 1 Example of non-rectangular domain Ω (solid dots)
 0 0 0 0 0 0 −1 −1
inside a 3 × 3 grid. When invoking the continuous optimality 0 0 0 0 0 0 1 1
condition, the discrete approximations of the Laplacian and
−1 0 0 −1 0 0 0 0
 
of the divergence near the boundary involve several points
inside ∂Ω (circles) for which no data is available. First-order  0 −1 0 0 −1 0 0 0 
approximation of the natural boundary condition (9) is thus
 
 0 0 −1 0 0 −1 0 0 
required. Relying only on discrete optimization simplifies a  
lot the boundary handling.
1 1 0 0 0 0 0 −1 0 
Dv =   (35)
2  0 1 0 0 0 0 0 −1 
0 0 1 0 0 1 0 0
+/−  
The sets Ωu/v all contain five pixels: 0 0 0 1 0 0 1 0
0 0 0 0 1 0 1 1
Ωu+ = {(1, 1) , (2, 1) , (1, 2) , (2, 2) , (1, 3)} (28)
Ωu− = {(2, 1) , (3, 1) , (2, 2) , (3, 2) , (2, 3)} (29)
Let us now show how these matrices relate to the
Ωv+ = {(1, 1) , (2, 1) , (3, 1) , (1, 2) , (2, 2)} (30) discretization of the continuous optimality condition (8).
Ωv− = {(1, 2) , (2, 2) , (3, 2) , (1, 3) , (2, 3)} (31) Using second-order central finite differences approxi-
+/−
mations of the Laplacian (∆zu,v ≈ zu,v−1 + zu−1,v +
so that the differentiation matrices Du/v have five non- zu+1,v + zu,v+1 − 4zu,v ) and of the divergence operator
zero rows according to their definition (17). For in- (∇ · gu,v ≈ 21 (pu+1,v − pu−1,v ) + 12 (qu,v+1 − qu,v−1 )),
stance, the matrix associated with the forward finite we obtain:
differences operator ∂u+ reads:

−1 1 0 0 0 0 0 0
 [4zu,v −zu,v−1 −zu−1,v −zu+1,v −zu,v+1 ]+λu,v zu,v =
 0 −1 1 0 0 0 0 0 1 1 0

 0 0 0 0 0 0 0 0
 [pu−1,v − pu+1,v ]+ [qu,v−1 − qu,v+1 ]+λu,v zu,v (36)
  2 2
+
 0 0 0 −1 1 0 0 0
Du =    (32) The pixel (u, v) = (2, 2) is the only one whose four
 0 0 0 0 −1 1 0 0

 0 0 0 0 0 0 0 0
  neighbors are inside Ω. In that case, (36) becomes:
 0 0 0 0 0 0 −1 1
0 0 0 0 0 0 0 0 [4z2,2 − z2,1 − z1,2 − z3,2 − z2,3 ] + λ2,2 z2,2
| {z } | {z }
The negative Laplacian matrix L defined in (23) is =(L)5,· z =(Λ2 )5,· z

worth: 1 1 0
= [p1,2 − p3,2 ] + [q2,1 − q2,3 ] + λ2,2 z2,2 (37)

2 −1 0 −1 0 0 0 0
 2
| {z } | 2 {z } | {z }
=(Λ2 )5,· z0
−1 3 −1 0 −1 0 0 0  =(Du )5,· p =(Dv )5,· q
 
 0 −1 2 0 0 −1 0 0 
 
−1 0 0 3 −1 0 −1 0  where we recognize the fifth equation of the discrete
L=   (33) optimality condition (26). This shows that, for pixels
 0 −1 0 −1 4 −1 0 −1

 0 0 −1 0 −1 2 0 0 
  having all four neighbors inside Ω, both the continuous
 0 0 0 −1 0 0 2 −1 and the discrete variational formulations yield the same
0 0 0 0 −1 0 −1 2 discretizations.
8 Yvain Quéau et al.

Now, let us consider a pixel near the boundary, for Using a similar rationale, we obtain equivalence of
instance pixel (1, 1). Using the same second-order dif- both formulations for the eight points inside Ω. Yet, let
ferences, (36) reads: us emphasize that discretizing the continuous optimal-
ity condition requires treating, on this example with a
[4z1,1 − z1,0 − z0,1 − z2,1 − z1,2 ] + λ1,1 z1,1 rather “simple” shape for Ω, not less than seven dif-
1 1 0 ferent cases (only pixels (3, 2) and (2, 3) are similar).
= [p0,1 − p2,1 ] + [q1,0 − q1,2 ] + λ1,1 z1,1 (38)
2 2 More general shapes bring out to play even more par-
which involves the values z1,0 and z0,1 of the depth map, ticular cases (points having only one neighbor inside Ω).
which we are not willing to estimate, and the values p0,1 Furthermore, boundary conditions must be invoked in
and q1,0 of the gradient field, which are not provided as order to approximate the depth values and the data out-
data. To eliminate these four values, we need to resort side Ω. On the other hand, the discrete functional pro-
to boundary conditions on z, p and q. The discretiza- vides exactly the same optimality condition, but with-
tions, using first order forward finite differences, of the out these drawbacks. The boundary conditions can be
natural boundary condition (9), at locations (1, 0) and viewed as implicitly enforced, hence PFreeB is satisfied.
(0, 1), read:
3.5 Empirical Evaluation
z1,1 − z1,0 = q1,0 (39)
z1,1 − z0,1 = p0,1 (40) We first consider the smooth surface from Figure 2,
whose normals are analytically known [26], and com-
hence the unknown depth values z1,0 and z0,1 can be pare three discrete least-squares methods which all sat-
eliminated from Equation (38): isfy PFast , PRobust and PFreeB : the DCT solution [53],
the Sylvester equations method [26], and the proposed
[2z1,1 − z2,1 − z1,2 ] + λ1,1 z1,1
one. As shown in Figures 2 and 3, our solution is slightly
1 1 0 more accurate. Indeed, the bias near the boundary in-
= [−p0,1 − p2,1 ] + [−q1,0 − q1,2 ] + λ1,1 z1,1 (41)
2 2 duced by the DCT method is corrected. On the other
Eventually, the unknown values p0,1 and q1,0 need to hand, we believe the reason why our method is more
be approximated. Since we have no information at all accurate than that from [26] is because we use a combi-
about the values of g outside Ω, we use homogeneous nation of forward and backward finite differences, while
Neumann boundary conditions13 : [26] relies on central differences. Indeed, when using
central differences to discretize the gradient, the second-
∇p · η = 0 over ∂Ω (42) order operator (Laplacian) appearing in the Sylvester
∇q · η = 0 over ∂Ω (43) equations from [26] involves none of the direct neigh-
bors, which may be non-robust for noisy data (see, for
Discretizing these boundary conditions using first order instance, Appendix 3 in [4]). For instance, let us con-
forward finite differences, we obtain: sider a 1D domain Ω with 7 pixels. Then, the following
differentiation matrix is advocated in [26]:
p0,1 = p1,1 (44)
−3 4 −1 0 0 0 0
 
q1,0 = q1,1 (45) −1 0 1 0 0 0 0
 
 0 −1 0 1 0 0 0
Using these identifications, the discretized optimality 1  
condition (41) is given by: Du =  0 0 −1 0 1 0 0 (47)
2 
 0 0 0 −1 0 1 0

[2z1,1 − z2,1 − z1,2 ] + λ1,1 z1,1  0 0 0 0 −1 0 1
| {z
=(L)1,· z
} | {z }
=(Λ2 )1,· z
0 0 0 0 1 −4 3
1 1 0
The optimality condition (Sylvester equation) in [26]
= [−p1,1 − p2,1 ] + [−q1,1 − q1,2 ] + λ1,1 z1,1 (46) involves the following second-order operator Du > Du :
2
| {z } | 2 {z } | {z }
=(Λ2 )1,· z0
10 −12 2 0 0 0
 
=(Du )1,· p =(Dv )1,· q 0
−12 17 −4 −1 0 0 0 
which is exactly the first equation of the discrete opti- 
 2 −4 3 0 −1 0 0 

mality condition (26). 1
Du > Du = 
 
0 -1 0 2 0 -1 0  (48)
4 
13
This assumption is weaker than the homogeneous Neu-
 0
 0 −1 0 3 −4 2 

mann boundary condition ∇z · η = 0 used by Agrawal et al.
 0 0 0 −1 −4 17 −12
in [1]. 0 0 0 0 2 −12 10
Variational Methods for Normal Integration 9

20 20

−20 −20

20 60 20 60
40 40 40 40
20 20
x 60 x 60
y y
Ground-truth Simchony et al. [53]

20 20

−20 −20

20 60 20 60
40 40 40 40
20 20
x 60 x 60
y y
Harker and O’Leary [26] Proposed

Fig. 2 Qualitative evaluation of the PRobust property. An additive, zero-mean, Gaussian noise with standard deviation
0.1kgk∞ was added to the (analytically known) gradient of the ground-truth surface, before integrating this gradient by
three least-squares methods. Ours qualitatively provides better results than the Sylvester equations method from Harker and
O’Leary [26]. It seems to provide similar robustness as the DCT solution from Simchony et al. [53], but the quantitative
evaluation from Figure 3 shows that our method is actually more accurate.

The bolded values of this matrix indicate that com- In addition, as predicted by the complexity analy-
putation of the second-order derivatives for the fourth sis in Subsection 3.2, our solution relying on precon-
pixel does not involve the third and fifth pixels. On ditioned conjugate gradient iterations has an asymp-
the other hand, with the proposed operator defined totic complexity (O(5n log(n) log(1/))) which is in-
in Equation (24), the second-order operator always in- between that of the Sylvester equations approach [26]
volves the “correct” neighborhood: (O(n1.5 )) and of DCT [53] (O(n log(n))). The CPU
times of our method and of the DCT solution, measured
1 −1 0 0 0 0 0
 
−1 2 −1 0 0 0 0  using Matlab codes running on a recent i7 processor, ac-

 0 −1 2 −1 0 0 0 
 tually seem proportional: according to this complexity
analysis, we guess the proportionality factor is around
Du > Du = 
 
 0 0 -1 2 -1 0 0  (49)

 0 0 0 −1 2 −1 0  5 log(1/). Indeed, with  = 10−4 , which is the value

 0 0 0 0 −1 2 −1
 we used in our experiments, 5 log(1/) ≈ 46, which is
consistent with the second graph in Figure 3.
0 0 0 0 0 −1 1
10 Yvain Quéau et al.

Simchony et al.
Harker and O’Leary
1.5 Proposed
RMSE (px)

RMSE = 0.11
0.5

0
0 0.05 0.1 0.15 0.2 0.25
σ
Simchony et al.
2
10 Harker and O’Leary
Proposed
CPU (s)

0 RMSE = 4.66
10

Fig. 4 3D-reconstruction of surface Svase (see Figure 3

in [48]) from its (analytically known) normals, using the
−2
10 proposed discrete least-squares method. Top: when Ω is re-
stricted to the image of the vase. Bottom: when Ω is the
whole rectangular grid. Quadratic integration smooths the
24

6
8

depth discontinuities and produces Gibbs phenomena near

04

09
12

25

51

10

4
8x

6x

2x

the kinks.
x

8x

6x
24
12

25

51

9
10

20

40

|Ω| 4 Piecewise Smooth Surfaces

We now tackle the problem of recovering a surface which

Fig. 3 Quantitative evaluation of the PRobust (top) and
is smooth only almost everywhere, i.e. everywhere ex-
PFast (bottom) properties. Top: RMSE between the depth
ground-truth and the ones reconstructed from noisy gradients cept on a “small” set where discontinuities and kinks
(adding a zero-mean Gaussian noise with standard deviation are allowed. Since all the methods discussed hereafter
σkgk∞ , for several values of σ). Bottom: Computation time rely on the same discretization as in Section 3, they in-
as a function of the size |Ω| of the reconstruction domain Ω.
herit its PFreeB and PNoRect properties, which will not
The method we put forward has a complexity which is inbe-
tween those of the methods of Simchony et al. [53] (based on be discussed in this section. Instead, we focus on the
DCT) and of Harker and O’Leary [26] (based on Sylvester PFast , PRobust , PNoPar , and of course PDisc properties.
equations), while being slightly more accurate than both of
them.

4.1 Recovering Discontinuities and Kinks

Besides its improved accuracy, the major advan- In order to clarify which variational formulations may
tage of our method over [26,53] is its ability to handle provide robustness to discontinuities, let us first con-
non-rectangular domains (PNoRect ). This makes possi- sider the 1D-example of Figure 5, with Dirichlet bound-
ble the 3D-reconstruction of piecewise-smooth surfaces, ary conditions. As illustrated in this example, least-
provided that a user segments the domain into pieces squares integration of a noisy normal field will provide
where z is smooth beforehand (see Figure 4). Yet, if a smooth surface. Replacing the least-squares estimator
the segmentation is not performed a priori , artifacts ΦL2 (s) = s2 by the sparsity one ΦL0 (s) = 1 − δ(s) will
are visible near the discontinuities, which get smoothed, minimize the cardinality of the difference between g and
and Gibbs phenomena appear near the continuous, yet ∇z, which provides a surface whose gradient is almost
non-differentiable kinks. We will discuss in the next sec- everywhere equal to g. As a consequence, robustness to
tion several strategies for removing such artifacts. noise is lost, yet discontinuities may be preserved.
Variational Methods for Normal Integration 11

10

Ground truth 9 ΦL0 (s)

Least-squares ΦL1 (s)
8
ΦL2 (s)
Sparsity
7 Φ1 (s)
Φ2 (s)
6

Φ(s)
5

4
Fig. 5 1D-illustration of integration of a noisy normal field
(arrows) over a regular grid (circles), in the presence of discon-
3
tinuities. The least-squares approach is robust to noise, but
2
smooths the discontinuities. The sparsity approach preserves
the discontinuities, but is not robust to noise. An ideal inte- 1
gration method would inherit robustness from least-squares,
and the ability to preserve discontinuities from sparsity. 0
−4 −2 0 2 4
s
These estimators can be interpreted as follows: least- Fig. 6 Graph of some robust estimators. The ability of ΦL2
to handle noise (small residuals) comes from its over-linear
squares assume that all residuals defined by k∇z(u, v)− behavior around zero, while that of ΦL0 to preserve disconti-
g(u, v)k are “low”, while sparsity assumes that most nuities (large residuals) is induced by its sub-linear behavior
of them are “zero”. The former is commonly used for in +∞. An estimator holding both these properties is neces-
“noise”, and the latter for “outliers”. In the case of nor- sarily non-convex (e.g., Φ1 and Φ2 , whose graphs are shown
with β = γ = 1), although ΦL1 may be an acceptable convex
mal integration, outliers may occur when: 1) ∇z(u, v) compromise.
exists but its estimate g(u, v) is not reliable; 2) ∇z(u, v)
is not defined because (u, v) lies within the vicinity of
a discontinuity or a kink. Considering that situation 1)
should rather be handled by robust estimation of the
gradient [31], we deal only with the second one, and use
the terminology “discontinuity” instead of “outlier”, al-
Another strategy consists in keeping least-squares as
though this also covers the concept of “kink”.
basis, but using it in a non-uniform manner. The sim-
We are looking for an estimator which combines plest way would be to remove the discontinuity points
the robustness of least-squares to noise, and that of from the integration domain Ω, and then to apply our
sparsity to discontinuities. These abilities are actually quadratic method from the previous section, since it
due to their asymptotic behaviors. Robustness of least- is able to manage non-rectangular domains. Yet, this
squares to noise comes from the quadratic behavior would require detecting the discontinuities beforehand,
around 0, which ensures that “low” residuals are consid- which might be tedious. It is actually more convenient
ered as “good” estimates, while this quadratic behav- to introduce weights in the least-squares functionals,
ior becomes problematic in ±∞: discontinuities yield which are inversely proportional to the probability of ly-
“high” residuals, which are over-penalized. The spar- ing on a discontinuity [47,50]. We discuss this weighted
sity estimator has the opposite behavior: treating the least-squares approach in Subsection 4.4, where a statis-
high residuals (discontinuities) exactly as the low ones tical interpretation of the Perona and Malik’s anisotropic
ensures that discontinuities are not over-penalized, yet diffusion model [44] is also exhibited. Eventually, an ex-
low residuals (noise) are. A good estimator would thus treme case of weighted least-squares consists in using bi-
be quadratic around zero, but sub-linear around ±∞. nary weights, where the weights indicate the presence
Obviously, only non-convex estimators hold both these of discontinuities. This is closely related to Mumford
properties. We will discuss several choices “inbetween” and Shah’s segmentation method [37], which simulta-
the quadratic estimator ΦL2 and the sparsity one ΦL0 neously estimates the discontinuity set and the surface.
(see Figure 6): the convex compromise ΦL1 (s) = |s| is We show in Subsection 4.5 that this approach is the
studied in Subsection 4.2, and the non-convex estima- one which is actually the most adapted to the problem
2
tors Φ1 (s) = log(s2 + β 2 ) and Φ2 (s) = s2s+γ 2 , where β of integrating a noisy normal field in the presence of
and γ are hyper-parameters, in Subsection 4.3. discontinuities.
12 Yvain Quéau et al.

4.2 Total Variation-like Integration The discrete functional to minimize is thus given by:

The problem of handling outliers in a noisy normal field 1 X X q 2  2

ETV (z) = ∂u+ zu,v −pu,v + ∂v+ zu,v −qu,v
has been tackled by Du, Robles-Kelly and Lu, who com- 4
(u,v)∈Ω ++
pare in [18] the performances of several M-estimators. X X q 2  2
They conclude that regularizers based on the L1 norm + ∂u+ zu,v −pu,v + ∂v− zu,v −qu,v
are the most effective ones. We provide in this sub- (u,v)∈Ω +−
section several numerical considerations regarding the X X q 2  2
+ ∂u− zu,v −pu,v + ∂v+ zu,v −qu,v
discretization of the L1 fidelity term:
(u,v)∈Ω −+
!
ZZ X X q 2  − 2
−
FL1 (z) = k∇z(u, v) − g(u, v)k1 du dv + ∂u zu,v −pu,v + ∂v zu,v −qu,v
(u,v)∈Ω −−
(u,v)∈Ω
XX  0
2
+ λu,v zu,v − zu,v (53)
ZZ n
= |∂u z(u, v) − p(u, v)|
(u,v)∈Ω
(u,v)∈Ω
o Minimizing (53) comes down to solving the following
+ |∂v z(u, v) − q(u, v)| du dv (50) constrained optimization problem:
1 XX XX
min krU V
u,v k
When p(u, v) ≡ 0 and q(u, v) ≡ 0, (50) is the so- z,{rU V } 4
(U,V )∈{+,−}2 (u,v)∈Ω U V
called “anisotropic total variation” (anisotropic TV) XX 2
0

regularizer, which tends to favor piecewise-constant so- + λu,v zu,v − zu,v
lutions while allowing discontinuity jumps. Considering (u,v)∈Ω
the discontinuities and kinks as the equivalent of edges s.t. rU V
u,v = ∇U V zu,v − gu,v (54)
in image restoration, it seems natural to believe that
the fidelity term (50) may be useful for discontinuity- where we denote ∇U V = [∂uU , ∂vV ]> , (U, V ) ∈ {+, −}2 ,
preserving integration. the discrete approximation of the gradient correspond-
This fidelity term is not only convex, but also de- ing to domain Ω U V .
couples the two directions u and v, which allows fast
ADMM-based (Bregman iterations) numerical schemes Numerical Solution. We solve the constrained optimiza-
involving shrinkages [24, 47]. On the other hand, it is tion problem (54) by the augmented Lagrangian method,
not so natural to use such a decoupling: if the value through an ADMM algorithm [21] (see [9] for a recent
of p is not reliable at some point (u, v), usually that overview of such algorithms). This algorithm reads:
of q is not reliable either. Hence, it may be wortwhile α XX XX
z(k+1) = argmin ∇U Vzu,v
to use instead a regularizer adapted from the “isotropic z∈R |Ω| 8
(U,V )∈{+,−}2 (u,v)∈Ω U V
TV”. This leads us to adapt the well-known model from   2
V (k) U V (k)
Rudin, Osher and Fatemi [49] to the integration prob- − gu,v +rU u,v −bu,v
lem: XX 2
0

+ λu,v zu,v − zu,v (55)
ZZ
(u,v)∈Ω
ETV (z) = k∇z(u, v) − g(u, v)k 2
V (k+1) α 
V (k)

(u,v)∈Ω rU
u,v = argmin r− ∇U Vzu,v
(k+1)
−gu,v +bU
u,v
r∈R2 8
2
+ λ(u, v) z(u, v) − z 0 (u, v) du dv (51)

+ krk (56)
V (k+1) V (k) V (k+1)
bU
u,v = bU
u,v +∇ UV (k+1)
zu,v − gu,v − rU
u,v
Discretization. Since the term k∇z(u, v) − g(u, v)k can (57)
be interpreted in different manners, depending on the
where the bU V are the scaled dual variables, and α > 0
neighborhood of (u, v), we need to discretize it appro-
corresponds to a descent stepsize, which is supposed
priately. Let us consider all four possible first-order dis-
to be fixed beforehand. Note that the choice of this
cretizations of the gradient ∇z, associated to the four
parameter influences only the convergence rate, not the
following sets of pixels:
actual minimizer. In our experiments, we used α = 1.
The z-update (55) is a linear least-squares problem
Ω U V = ΩuU ∩ ΩvV , (U, V ) ∈ {+, −}2 (52) simimilar to the one which was tackled in Section 3.
Variational Methods for Normal Integration 13

σ = 0% - RMSE = 4.52 σ = 0.5% - RMSE = 4.62 σ = 1% - RMSE = 4.79

Fig. 7 Depth estimated after 1000 iterations of the TV-like approach, in the presence of additive, zero-mean, Gaussian noise
with standard deviation equal to σkgk∞ . The indicated RMSE is computed on the whole domain. In the absence of noise,
both discontinuities and kinks are restored, although staircasing artifacts appear. In the presence of noise, the discontinuities
are smoothed. Yet, the 3D-reconstruction near the kinks is still more satisfactory than the least-squares one: Gibbs phenomena
are not visible, unlike in the second row of Figure 4.

Its solution z(k+1) is the solution of the following SDD can be improved rather simply [23]. We cannot con-
linear system: sider that PFast is satisfied since, in comparison with the
(k)
quadratic method from Section 3, yet the TV approach
ATV z(k+1) = bTV (58) is “reasonably” fast. Possibly faster algorithms could
be employed, as for instance the FISTA algorithm from
with :
Beck and Teboulle [7], or primal-dual algorithms [13],
α X Xh U> U >
i
but we leave such improvements as future work.
ATV = Du Du + DVv DVv +Λ2 (59)
8 On the other hand, according to the results from
(U,V )∈{+,−}2
(k) α X X h U > U V (k) > (k)
i Figure 7, discontinuities are recovered in the absence of
bTV = Du p + DVv qU V +Λ2 z0 (60) noise, although staircasing artifacts appear (such arti-
8 2
(U,V )∈{+,−}
facts are partly due to the non-differentiability of TV in
U/V zero [38]). Yet, the recovery of discontinuities is deceiv-
where the Du/v matrices are defined as in (17), the
ing when the noise level increases. On noisy datasets,
(k)
Λ matrix as in (20), and where we denote pU V and the only advantage of this approach over least-squares
(k) (k) (k)
qU V the components of g + rU V − bU V . is thus that it removes the Gibbs phenomena around
The solution of System (58) can be approximated the kinks i.e., where the surface is continuous, but non-
by conjugate gradient iterations, choosing at each iter- differentiable (e.g., the sides of the vase).
ation the previous estimate z(k) as initial guess (setting Because of the staircasing artifacts and of the lack of
z(0) , for instance, as the least-squares solution from Sec- robustness to noise, we cannot find this first approach
tion 3). In addition, the matrix ATV is always the same: satisfactory. Yet, since turning the quadratic functional
this allows computing the preconditioner only once. into a non-quadratic one seems to have positive influ-
Eventually, the r-updates (56), (u, v) ∈ Ω, are basis ence on discontinuities recovery, we believe that explor-
pursuit problems [17], which admit the following closed- ing non-quadratic models is a promising route. Stair-
form solution (generalized shrinkage): casing artifacts could probably be reduced by replacing
total variation by total generalized variation [10], but
 U V (k+1)
su,v we rather consider now non-convex models.

(k+1) V (k+1) 4
rU V
u,v = max ksU
u,v k − , 0 (61)
α ksU V (k+1) k
u,v
4.3 Non-convex Regularization
with:
(k+1) (k) Let us now consider non-convex estimators Φ in the
sU V
u,v = ∇U V zu,v
(k+1)
− gu,v + bU V
u,v (62)
fidelity term (5), which are often referred to as “Φ-
Discussion. This TV-like approach has two main ad- functions” [4]. As discussed in Subsection 4.1, the choice
vantages: apart from the stepsize α which controls the of a specific Φ-function should be made according to
speed of convergence, it does not depend on the choice several principles:
of a parameter, and it is convex. The initialization has • Φ should have a quadratic behavior around zero,
influence only on the speed of convergence, and not on in order to ensure that the integration is guided by
the actual minimizer: convergence towards the global the “good” data. The typical choice ensuring this
minimum is guaranteed [51]. It can be shown that the property is ΦL2 (s) = s2 , which was discussed in
convergence rate of this scheme is ergodic, and this rate Section 3;
14 Yvain Quéau et al.

• Φ should have a sublinear behavior at infinity, so Introducing the notations:

that outliers do not have a predominant influence, 1 XX XX
Φ k∇U Vzu,v − gu,v k (66)


and also to preserve discontinuities and kinks. The f (z)=
4
typical choice is the sparsity estimator ΦL0 (s) = 0 2 UV
(U,V )∈{+,−} (u,v)∈Ω
if s = 0 and ΦL0 (s) = 1 otherwise; g(z) = kΛ(z − z0 )k2 (67)
• Φ should ideally be a convex function.
the discrete functional (65) is rewritten:
Obviously, it is not possible to simultaneously sat-
isfy these three properties. The TV-like fidelity term EΦ (z) = f (z) + g(z) (68)
introduced in Subsection 4.2 is a sort of “compromise”:
where f is smooth, but non-convex, and g is convex
it is the only convex function being (over-) linear in 0
(and smooth, although non-smooth functions g could
and (sub-) linear in ±∞. Although it does not depend
be handled).
on the choice of any hyper-parameter, we saw that it
has the drawback of yielding the so-called “staircase
Numerical Solution. The problem of minimizing a dis-
effect”, and that discontinuities were not recovered so
crete energy like (68), yielded by the sum of a con-
well in the presence of noise. If we accept to lose the
vex term g and a non-convex, yet smooth term f , can
convexity of Φ, we can actually design estimators which
be handled by forward-backward splitting. We use the
better fit both other properties. Although there may
“iPiano” iterative algorithm by Ochs et al. [41], which
then be several minimizers, such non-convex estima-
reads:
tors were recently shown to be very effective for image   
restoration [36]. −1
z(k+1)= (I+α1 ∂g) z(k)−α1∇f (z(k) )+α2 z(k)−z(k−1)
We will consider two classical Φ-functions, whose
graphs are plotted in Figure 6: (69)

 2s where α1 and α2 are suitable descent stepsizes (in our

Φ01 (s) = 2

Φ1 (s) = log(s2 + β 2 )  implementation, α2 is fixed to 0.8, and α1 is chosen by
s + β2
s2 ⇒ 2 (63) the “lazy backtracking” procedure described in [41]),
Φ2 (s) = Φ02 (s) = 2 γ s
 −1
(I + α1 ∂g) is the proximal operator of g, and ∇f (z(k) )
s2 + γ 2 (s + γ 2 )2
2
is the gradient of f evaluated at current estimate z(k) .
Let us remark that these estimators were initially We detail hereafter how to evaluate the proximal oper-
introduced in [19] in this context, and that other non- ator of g and the gradient of f .
convex estimators can be considered, based for instance The proximal operator of g writes, using (67):
on Lp norms, with 0 < p < 1 [5]. kx − xbk
−1
Let us now show how to numerically minimize the (I + α1 ∂g) (b
x) = argmin + α1 g(x) (70)
x∈R|Ω| 2
resulting functionals:
−1
= I + 2α1 Λ2 b + 2α1 Λz0 (71)


ZZ x
EΦ (z) = Φ (k∇z(u, v) − g(u, v)k) where the inversion is easy to compute, since the matrix
(u,v)∈Ω involved is diagonal.
 0
2 In order to obtain a closed-form expression of the
+ λ(u, v) z(u, v) − z (u, v) du dv (64)
gradient of f defined in (66), let us rewrite this function
in the following manner:
Discretization. We consider the same discretization strat-
egy as in Subsection 4.2, aiming at minimizing the dis- 1 XX XX
Φ kDU V


f (z) = u,v z − gu,v k (72)
crete functional: 4 2 UV
(U,V )∈{+,−} (u,v)∈Ω

1 XX XX  
where DU V
EΦ (z) = Φ ∇U Vzu,v − gu,v u,v is a 2 × |Ω| finite differences matrix used
4 for approximating the gradient at location (u, v), using
(U,V )∈{+,−}2 (u,v)∈Ω U V
XX  0
2 the finite differences operator ∇U V , {U, V } ∈ {+, −}2 :
+ λu,v zu,v − zu,v (65)
(u,v)∈Ω "
#
DUu m−1 (u,v),·
which resembles the TV functional defined in (53), and DU V
u,v =

(73)
DVv m−1 (u,v),·
where ∇U V represents the finite differences approxima-
tion of the gradient used over the domain Ω U V , with where we recall that the mapping m associates linear
{U, V } ∈ {+, −}2 . indices with pixel coordinates (see Equation (18)).
Variational Methods for Normal Integration 15

β = 0.1 - RMSE = 4.60 β = 0.5 - RMSE = 4.42 β = 1 - RMSE = 5.08

γ = 0.5 - RMSE = 4.51 γ = 1 - RMSE = 4.44 γ = 5 - RMSE = 4.67

Fig. 8 Non-convex 3D-reconstructions of surface Svase , using Φ1 (top) or Φ2 (bottom). An additive, zero-mean, Gaussian
noise with standard deviation σkgk∞ , σ = 1%, was added to the gradient field. The non-convex approaches depend on the
tuning of a parameter (β or γ), but they are able to reconstruct the discontinuities in the presence of noise, unlike the TV
approach. Staircasing artifacts indicate the presence of local minima (we used as initial guess z (0) the least-squares solution).

The gradient of f is thus given by: of a “Canadian tent”-like surface, with additive, zero-
( mean, Gaussian noise (σ = 10%), is presented. When
1 XX XX
V> using the least-squares solution as initial guess z (0) , the
DU DU V


∇f (z) = u,v u,v z − gu,v
4 3D-reconstruction is very close to the genuine surface.
(U,V )∈{+,−}2 (u,v)∈Ω U V

) Yet, when using the trivial initialization z (0) ≡ 0, we
Φ0 kDU V
u,v z − gu,v k obtain a surface whose slopes are “almost everywhere”
× (74)
kDU V
u,v z − gu,v k equal to the real ones, but unexpected discontinuity
jumps appear. Since only the initialization differs in
Given the choices (63) for the Φ-functions, this can be these experiments, this clearly shows that the artifacts
further simplified: indicate the presence of local minima.

X X DU V>
DU V

Although local minima can sometimes be avoided
1 XX u,v u,v z − gu,v
∇f1 (z) = by using the least-squares solution as initial guess (e.g.,
2 kDU V 2
u,v z − gu,v k + β
2
(U,V )∈{+,−}2 (u,v)∈Ω
UV Figure 9), this is not always the case (e.g., Figure 8).
(75) Hence, the non-convex estimators perform overall bet-
X X γ 2 DU V>

ter than the TV-like approach, but they are still not
1 XX u,v DU V
u,v z − gu,v
∇f2 (z) =
optimal. We now follow other routes, which use least-
2 UV 2 2 2
(U,V )∈{+,−}2 (u,v)∈Ω U V kDu,v z − gu,v k + γ squares as basis estimator, yet in a non-uniform man-
(76) ner, in order to allow discontinuities.

Discussion. Contrarily to the TV-like approach (see

Subsection 4.2), the non-convex estimators require set-
ting one hyper-parameter (β or γ). As shown in Fig-
ure 8, the choice of this parameter is crucial: when it is
too high, discontinuities are smoothed, while setting a 4.4 Integration by Anisotropic Diffusion
too low value leads to strong staircasing artifacts. Inbe-
tween, the values β = 0.5 and γ = 1 seem to preserve Both previous methods (total variation and non-convex
discontinuities, even in the presence of noise (which was estimators) replace the least-squares estimator by an-
not the case using the TV-like approach). other one, assumed to be robust to discontinuities. Yet,
Yet, staircasing artifacts are still present. Despite it is possible to proceed differently: the 1D-graph in
their non-convexity, the new estimators Φ1 and Φ2 are Figure 5 shows that most of data are corrupted only by
differentiable, hence these artifacts do not come from noise, and that the discontinuity set is “small”. Hence,
a lack of differentiability, as this was the case for TV. applying least-squares everywhere except on this set
They rather indicate the presence of local minima. This should provide an optimal 3D-reconstruction. To achieve
is illustrated in Figure 9, where the 3D-reconstruction this, a first possibility is to consider weighted least-
16 Yvain Quéau et al.

the continuous optimality condition associated to (77)

is related to their “anisotropic diffusion model” 14 . Such
tensor fields W : Ω → R2×2 are called “diffusion ten-
sors”: we refer the reader to [54] for a complete overview.
The use of diffusion tensors for the integration prob-
lem is not new [47], but we provide hereafter additional
comments on the statistical interpretation of such ten-
sors. Interestingly, the diffusion tensor (78) also appears
Ground-truth
when making different assumptions on the noise model
than those we considered so far. Up to now, we assumed
that the input gradient field g was equal to the gradient
∇z of the depth map z, up to an additive,
 zero-mean,
 2 
> σ 0
Gaussian noise: g = ∇z + ,  ∼ N [0, 0] , .
0 σ2
This hypothesis may not always be realistic. For in-
stance, in 3D-reconstruction scenarii such as photomet-
ric stereo [55], one estimates the normal field n : Ω →
z (0) = least-squares solution - RMSE = 0.78 R3 pixelwise, rather than the gradient g : Ω → R2 ,
from a set of images. Hence, the Gaussian assumption
should rather be made on these images. In this case,
and provided that a maximum-likelihood for the nor-
mals is used, it may be assumed that the estimated nor-
mal field is the genuine one, up to an additive Gaussian
noise. Yet, this does not imply that the noise in the
gradient field g is Gaussian-distributed. Let us clarify
this point.
z (0) ≡ 0 - RMSE = 13.16 Assuming orthographic projection, the relationship
>
between n = [n1 , n2 , n3 ] and ∇z is written, in every
Fig. 9 3D-reconstruction of a “Canadian tent”-like surface point (u, v) where the depth map z is differentiable:
from its noisy gradient (σ = 1%), by the non-convex inte-
grator Φ1 (β = 0.5, 12000 iterations), using two different
1 >
initializations. The objective function being non-convex, the −∇z(u, v)> , 1

n(u, v) = p (79)
iterative scheme may converge towards a local minimum. k∇z(u, v)k2 + 1

squares: which implies that [− nn31 , − nn23 ]> = [∂u z, ∂v z]> = ∇z. If
>
ZZ
2
we denote n = [n1 , n2 , n3 ] the estimated normal field,
min kW(u, v) [∇z(u, v) − g(u, v)]k it follows from (79) that [− nn13 , − nn23 ]> = [p, q]> = g.
z
(u,v)∈Ω Let us assume that n and n differ according to an
2 additive, zero-mean, Gaussian noise:
+ λ(u, v) z(u, v) − z 0 (u, v) du dv

(77)
where W is a Ω → R2×2 tensor field, acting as a weight n(u, v) = n(u, v) + (u, v) (80)
map designed to reduce the influence of discontinuity
points. The weights can be computed beforehand ac-
where :
cording to the integrability of g [47], or by convolution
of the components of g by a Gaussian kernel [1]. Yet,   2 
σ 0 0
such approaches are of limited interest when g contains
(u, v) ∼ N [0, 0, 0]> ,  0 σ 2 0  (81)
noise. In this case, the weights should rather be set as
0 0 σ2
a function inversely proportional to k∇z(u, v)k, e.g.:
1 Since n3 is unlikely to take negative values (this
W(u, v) = r 2 I2 (78)
k∇z(u,v)k would mean that the estimated surface is not oriented
µ +1
14
Although (78) actually yields an isotropic diffusion model,
with µ a user-defined hyper-parameter. The latter ten- since it “utilizes a scalar-valued diffusivity and not a diffusion
sor is the one proposed by Perona and Malik in [44]: tensor” [54].
Variational Methods for Normal Integration 17

towards the camera), the following Geary-Hinkley trans- Discretization. Using the same discretization strategy
forms: as in Subsections 4.2 and 4.3 leads us to the following
  discrete functional:
n3 nn13 − n1 (
t1 = s  (82) 1 XX
2
EPM (z) = AU V (z) DUu z−p
 2 
n 4
σ2 n3
1
+1 (U,V )∈{+,−}2
)
 
2
n2 UV V
n3 n3 − n2 + B (z) Dv z−q
t2 = s   (83)
2
n2 2
+ Λ z − z0


σ2 n3 +1 (89)

where the AU V (z) and BU V (z) are |Ω| × |Ω| diagonal

both follow standard Gaussian distribution N (0, 1) [27].
matrices containing the following values:
After some algebra, this can be rewritten as:
1
1 aU V
u,v = q (90)
[∂u z − p] ∼ N (0, 1) (84)
q 2 2
pu,v
2 Uz
(∂u V
u,v ) +(∂v zu,v )
1+ +1
p p
σ 1+p 2 k∇zk2 + 1 ν µ2
1 1
p p [∂v z − q] ∼ N (0, 1) (85) bU V
u,v = q q (91)
σ 1 + q 2 k∇zk2 + 1 qu,v 2

(∂ zu,v )2 +(∂ V zu,v )2
U
1+ ν
u
µ2
v
+1
This rationale suggests the use of the following fi-
with (U, V ) ∈ {+, −}2 .
delity term:
ZZ
FPM (z) =
2
kW(u, v) [∇z(u, v) − g(u, v)]k du dv Numerical Solution. Since the coefficients aU V UV
u,v and bu,v
depend in a nonlinear way on the unknown values zu,v ,
(u,v)∈Ω
it is difficult to derive a closed-form expression for the
(86)
minimizer of (89). To deal with this issue, we use the
where W(u, v) is the following 2×2 anisotropic diffusion following fixed point scheme, which iteratively updates
tensor field: the anisotropic diffusion tensors and the z-values:
" #
√ 1
(
1 0 1
2
1+p(u,v)2
XX
W(u, v) = p 0 √ 1 (87) z(k+1) = argmin AU V (z(k) ) DU u z−p
k∇z(u, v)k2 +1 1+q(u,v)2 z∈R|Ω| 4(U,V )∈{+,−}2
)
Unfortunately, we experimentally found with the
2
UV (k) V
choice (87) for the diffusion tensor field, discontinuities + B (z ) Dv z−q
were not always recovered. Instead, following the pio-
2
+ Λ z − z0


neering ideas from Perona and Malik [44], we introduce (92)
two parameters µ and ν to control the respective in-
fluences of the terms depending on the gradient of the Now that the diffusion tensor coefficients are fixed, each
unknown k∇zk and on the input gradient (p, q). The optimization problem (92) is reduced to a simple linear
new tensor field is then given by: least-squares problem. In our implementation, we solve
  the corresponding optimality condition using Cholesky
1
0 factorization, which we experimentally found to provide
1
q
p(u,v) 2
1+( )
W(u, v) = r 2

0
ν
q 1
 (88) more stable results than conjugate gradient iterations.
q(u,v) 2
k∇z(u,v)k 1+( )
µ +1 ν

Discussion. We first experimentally verify that the pro-

Replacing the matrix in (88) by I2 yields exactly the posed anisotropic diffusion approach is indeed a statis-
Perona-Malik diffusion tensor (78), which reduces the tically meaningful approach in the context of photomet-
influence of the fidelity term on locations (u, v) where ric stereo. As stated in [39], “in previous work on pho-
k∇z(u, v)k increases, which are likely to indicate dis- tometric stereo, noise is [wrongly] added to the gradi-
continuities. Yet, our diffusion tensor (88) also reduces ent of the height function rather than camera images”.
the influence of points where p or q is high, which are Hence, we consider the images from the “Cat” dataset
also likely to correspond to discontinuities. In our ex- presented in [52], and add a zero-mean, Gaussian noise
periments, we found that ν = 10 could always be used, with standard deviation σkIk∞ , σ = 5%, to the im-
yet the choice of µ has more influence on the actual ages, where kIk∞ is the maximum graylevel value. The
results. normals were computed by photometric stereo [55] over
18 Yvain Quéau et al.

45
the part representing the cat. Then, since only the nor-
Simchony et al.
mals ground-truth is provided in [52], and not the depth 40
Harker and O’Leary
ground-truth, we a posteriori computed the final nor- 35
Proposed (least-squares)

MAE (degrees)
Proposed (anisotropic diffusion)
mal maps by central finite differences. This allows us to 30
calculate the angular error, in degrees, between the real
25
surface and the reconstructed one. The mean angular
error (MAE) can eventually be computed over the set 20

of pixels for which central finite differences make sense 15

(boundary and background points are excluded). 10

0 0.05 0.1 0.15 0.2
σ

Fig. 11 Mean angular error (in degrees) as a function of

the standard deviation σkIk∞ of the noise which was added
to the photometric stereo images. The anisotropic diffusion
approach always outperforms least-squares. For the meth-
ods [26, 53], the gradient field was filled with zeros outside
the reconstruction domain, which adds even more bias.

Although the parameter-free diffusion tensor (87)

seems able to recover discontinuities, this is not always
the case. For instance, we did not succeed in recov-
ering the discontinuities of the surface Svase . For this
dataset, we had to use the tensor (88). The results from
Least-squares Anisotropic diffusion Figure 12 show that with an appropriate tuning of µ,
discontinuities are recovered and Gibbs phenomena are
removed, without staircasing artifact. Yet, as in the ex-
periment of Figure 10, the discontinuities are not very
sharp. Such artifacts were also observed by Badri et
al. [5], when experimenting with the anisotropic diffu-
sion tensor from Agrawal et al. [1]. Sharper discontinu-
ities could be recovered by using binary weights: this is
the spirit of the Mumford-Shah segmentation method,
which we explore in the next subsection.

(MAE = 9.29 degrees) (MAE = 8.43 degrees) 4.5 Adaptation of the Mumford and Shah Functional

Fig. 10 Top row: three out of the 96 input images used for Let z 0 : Ω → R be a noisy image to restore. In order to
estimating the normals by photometric stereo [55]. Middle estimate a denoised image z while perserving the dis-
row, left: 3D-reconstruction by least-squares integration of
the normals (see Section 3). Bottom row, left: angular error continuities of the original image, Mumford and Shah
map (blue is 0 degree, red is 60 degrees). The estimation is suggested in [37] to minimize a quadratic functional
biased around the occluded areas. Middle and bottom rows, only over a subset Ω\K of Ω, while automatically esti-
right: same, using anisotropic diffusion integration with the mating the discontinuity set K according to some prior.
tensor field defined in (87). The errors remain confined in the
occluded parts, and do not propagate over the discontinuities. A reasonable prior is that the length of K is “small”,
which leads to the following optimization problem:
ZZ Z
min µ k∇z(u, v)k2 du dv + dσ
Figure 10 shows that the 3D-reconstruction obtained z,K K
(u,v)∈Ω\K
by anisotropic diffusion outperforms that obtained by ZZ
least-square: discontinuities are partially recovered, and 2
z(u, v) − z 0 (u, v) du dv

+λ (93)
robustness to noise is improved (see Figure 11). How-
(u,v)∈Ω\K
ever, although the diffusion tensor (87) does not require R
any parameter tuning, the restoration of discontinuities where λ and µ are positive constants, and K dσ is the
is not as sharp as with the non-convex integrators, and length of the set K. See [4] for a detailed introduction
artifacts are visible along the discontinuities. to this model and its qualitative properties.
Variational Methods for Normal Integration 19

We modify the above models, so that they fit our

integration problem. Considering g as basis for least-
squares integration everywhere except on the disconti-
nuity set K, we obtain the following energy:
ZZ Z
EMS (z, K) = µ k∇z(u, v) − g(u, v)k2 du dv + dσ
K
(u,v)∈Ω\K
ZZ
2
λ(u, v) z(u, v)−z 0 (u, v) du dv (95)

µ = 0.02 - RMSE = 2.38 +
(u,v)∈Ω\K

for the Mumford-Shah functional, and the following

Ambrosio-Tortorelli approximation:
ZZ
EAT (z, w) = µ w(u, v)2 k∇z(u, v) − g(u, v)k2 du dv
(u,v)∈Ω
ZZ  
1
µ = 0.2 - RMSE = 2.19 +  k∇w(u, v)k2 + [w(u, v) − 1]2 du dv
4
(u,v)∈Ω
ZZ
2
λ(u, v) z(u, v) − z 0 (u, v) du dv

+ (96)
(u,v)∈Ω

where w : Ω → R is a smooth approximation of 1−χK .

Numerical Solution. We use the same strategy as in

µ = 2 - RMSE = 5.09
Section 3 for discretizing ∇z(u, v) inside Functional (96),
i.e. all the possible first-order discrete approximations
Fig. 12 Integration of the noisy gradient of Svase (σ = 1%)
of the differential operators are summed. Since disconti-
by anisotropic diffusion. As long as µ is small enough, discon- nuities are usually “thin” structures, it is possible that a
tinuities are recovered. Besides, no staircasing artifact is vis- forward discretization contains the discontinuity while
ible. Yet, the restored discontinuities are not perfectly sharp. a backward discretization does not. Hence, the defini-
tion of the weights w should be made accordingly to
+/−
Several approaches have been proposed to numeri- that of ∇z. Thus, we define four fields wu/v : Ω → R,
+/−
cally minimize the Mumford-Shah functional: finite dif- associated with the finite differences operators ∂u/v .
ferences scheme [12], piecewise constant approximation This leads to the following discrete analogue of Func-
[14], primal-dual algorithms [46], etc. Another approach tional (96):
consists in using elliptic functionals. An auxiliary func-
tion w : Ω → R is introduced. This function stands for EAT (z, wu+ , wu− , wv+ , wv− ) =
1 − χK , where χK is the characteristic function of the µ 2 2
Wu+ D+ + Wu− D−



set K. Ambrosio and Tortorelli have proposed in [2] to uz−p uz−p
2
consider the following optimization problem: !
2 2
Wv+ D+ Wv− D−



+ vz −q + vz −q
ZZ
min µ w(u, v)2 k∇z(u, v)k2 du dv
z,w
(u,v)∈Ω  + + 2 − 2 + 2 − 2

+ Du wu + D−u wu + D+
v wv + D−
v wv
ZZ 
1
 2
2 2
+  k∇w(u, v)k + [w(u, v)−1] du dv 1 + 2 2 2 2

4 + wu −1 + wu− −1 + wv+ − v1 + wv− −1
(u,v)∈Ω 8

2
+ Λ z − z0
ZZ
2 (97)
z(u, v) − z 0 (u, v) du dv

+λ (94)
+/−
(u,v)∈Ω where wu/v ∈ R|Ω| is a vector containing the values of
+/− +/− +/−
By using the theory of Γ -convergence, it is possible to the discretized field wu/v , and Wu/v = Diag(wu/v )
show that (94) is a way to solve (93) when  → 0. is the |Ω|×|Ω| diagonal matrix containing these values.
20 Yvain Quéau et al.

µ = 1 - RMSE = 4.94 µ = 45 - RMSE = 2.37 µ = 100 - RMSE = 4.14

Fig. 13 3D-reconstructions from the noisy gradient of Svase (σ = 1%), using the Mumford-Shah integrator. If µ is tuned
appropriately, sharp discontinuities can be restored, without staircasing artifacts.

We tackle the nonlinear problem (97) by an alter- The Mumford-Shah functional being non-convex, lo-
nating optimization scheme: cal minima may exist. Yet, as shown in Figure 14, the
choice of the initialization may not be as crucial as with
(k) (k) (k) (k)
z(k+1) = argmin EAT (z, wu+ ,wu− ,wv+ ,wv− ) (98) the non-convex estimators from Subsection 4.3. Indeed,
z∈R|Ω| the 3D-reconstruction of the “Canadian tent” surface is
(k+1) (k) (k) (k) similar using as initial guess the least-squares solution
wu+ = argmin EAT (z(k+1) , w,wu− ,wv+ ,wv− ) (99)
w∈R|Ω| or the trivial initialization z (0) ≡ 0.

and similar straightforward updates for the other in-

dicator functions. We can choose as initial guess, for
instance, the smooth solution from Section 3 for z(0) ,
(0) (0) (0) (0)
and wu+ = wu− = wv+ = wv− ≡ 1.
At each iteration (k), updating the surface and the
indicator functions requires solving a series of linear
least-squares problems. We achieve this by solving the
resulting linear systems (normal equations) by means
of the conjugate gradient algorithm. Contrarily to the
approaches that we presented so far, the matrices in- z (0) = least-squares solution - RMSE = 0.74
volved in these systems are modified at each iteration.
Hence, it is not possible to compute the preconditioner
beforehand. In our experiments, we did not consider
any preconditioning strategy at all. Thus, the proposed
scheme could obviously be accelerated.

Discussion. Let us now check experimentally, on the

same noisy gradient of surface Svase as in previous ex- z (0) ≡ 0 - RMSE = 1.84
periments, whether the Mumford-Shah integrator sat-
isfies the expected properties. In the experiment of Fig- Fig. 14 3D-reconstructions of the “Canadian tent” surface
ure 13, we performed 50 iterations of the proposed al- from its noisy gradient (σ = 1%), by the Mumford-Shah in-
ternating optimization scheme, with various choices for tegrator (µ = 20), using two different initializations. The ini-
the hyper-parameter µ. The  parameter was set to tialization matters, but not as much as with the non-convex
estimators from Subsection 4.3.
 = 0.1 (this parameter is not critical: it only has to
be “small enough”, in order for the Ambrosio-Tortorelli
approximation to converge towards the Mumford-Shah Hence, among all the variational integration meth-
functional). As it was already the case with other non- ods we have studied, the adaptation of the Mumford-
convex regularizers (see Subsection 4.3), a bad tuning Shah model is the approach which provides the most
of the parameter leads either to over-smoothing (low satisfactory 3D-reconstructions in the presence of sharp
values of µ) or to staircasing artifacts (high values of features: it is possible to recover discontinuities and
µ), which indicate the presence of local minima. Yet, kinks, even in the presence of noise, and with limited
by appropriately setting this parameter, we obtain a artifacts. Nevertheless, local minima may theoretically
3D-reconstruction which is very close to the genuine arise, as well as staircasing if the parameter µ is not
surface, without staircasing artifact. tuned appropriately.
Variational Methods for Normal Integration 21

Table 1 Main features of the five methods of integration proposed in this paper. The quadratic method has all desirable prop-
erties, except PDisc . The others lose PFast but hold PDisc . Sharpest features are recovered by using non-convex regularization
or the Mumford-Shah approach, yet staircasing artifacts and local minima may appear. In addition, all discontinuity-preserving
methods except TV require tuning at least one hyper-parameter. Yet, TV is not able to recover discontinuities in the presence
of noise. Overall, we recommend using: quadratic integration if speed is the most important issue; the Mumford-Shah approach
if recovering discontinuities is the most important issue; and anisotropic diffusion if discontinuities are present, but limited.

Method PFast PRobust PFreeB PDisc PNoRect PNoPar Local minima Staircasing
Quadratic +++ + + − + + No No
Total variation + + + + + + No Yes
Non-convex regularization − + + +++ + − Yes Yes
Anisotropic diffusion − + + ++ + −− No No
Mumford-Shah − + + +++ + −− Yes Yes

5 Conclusion and Perspectives integration. Quadratic (Tikhonov) smoothness terms

were suggested in [26]: to enforce surface smoothness
We proposed several new variational methods for solv- while preserving the discontinuities, we should rather
ing the normal integration problem. These methods were consider non-quadratic ones. In this view, higher-order
designed to satisfy the largest subset of properties that functionals (e.g., total generalized variation methods [10])
were identified in a companion survey paper [48] enti- may reduce not only these artifacts, but also staircas-
tled Normal Integration: A Survey. ing. Indeed, as shown in Figure 15, such artifacts may
We first detailed in Section 3 a least-squares solu- be visible when performing photometric stereo [55] with-
tion which is fast, robust and parameter-free, while as- out prior segmentation. Yet, this example also shows
suming neither a particular shape for the integration that the artifacts are visible only over the background,
domain nor a particular boundary condition. However, and do not seem to affect the relevant part.
discontinuities in the surface can be handled only if 3D-reconstruction is not the only application where
the integration domain is first segmented into pieces efficient tools for gradient field integration are required.
without discontinuities. Therefore, we discussed in Sec- Although the assumption on the noise distribution may
tion 4 several non-quadratic or non-convex variational differ from one application to another, PDE-based imag-
formulations aiming at appropriately handling discon- ing problems such as Laplace image compression [45] or
tinuities. As we have seen, the latter property can be Poisson image editing [43] also require an efficient inte-
satisfied only if (slow) iterative schemes are used and / grator. In this view, the ability of our methods to han-
or one critical parameter is tuned. Therefore, there is dle control points may be useful. We illustrate in Fig-
still room for improvement: a fast, parameter-free in- ure 16 an interesting application. From an RGB image
tegrator, able to handle discontinuities remains to be I, we selected the points where the norm of the gra-
proposed. dient of the luminance (in the CIE-LAB color space)
Table 1 summarizes the main features of the five was the highest (conserving only 10% of the points).
new integration methods proposed in this article. Con- Then, we created a gradient field g equal to zero ev-
trarily to Table 1 in [48], which recaps the features of erywhere, except on the control points, where it was
state-of-the-art methods, this time we use a more nu- set to the gradient of the color levels. The prior z 0 was
anced evaluation than binary features +/−. Among the set to a null scalar field, except on the control points
new methods, we believe that the least-squares method where we retained the original color data. Eventually,
discussed in Section 3 is the best if speed is the most λ is set to an arbitrary small value (λ = 10−9 ) ev-
important criterion, while the Mumford-Shah approach erywhere, except on the control points (λ = 10). The
discussed in Subsection 4.5 is the most appropriate one integration of each color channel gradient is performed
for recovering discontinuities and kinks. Inbetween, the independently, using the Mumford-Shah method to ex-
anisotropic diffusion approach from Subsection 4.4 rep- trapolate the data from the control points to the whole
resents a good compromise. grid. Using this approach, we obtain a nice piecewise-
Future research directions may include accelerating constant approximation of the image, in the spirit of
the numerical schemes and proving their convergence the “texture-flattening” application presented in [43].
when this is not trivial (e.g., for the non-convex inte- Besides, by selecting the control points in a more opti-
grators). We also believe that introducing additional mal way [8, 28], this approach could easily be extended
smoothness terms inside the functionals may be use- to image compression, reaching state-of-the-art lossy
ful for eliminating the artifacts in anisotropic diffusion compression rates. In fact, existing PDE-based meth-
22 Yvain Quéau et al.

(a) (b) (c) (d)

(e) (f)

Fig. 15 3D-reconstruction using photometric stereo. (a-c) All (real) input images. (d) 3D-reconstruction by least-squares on
the whole grid. (e) 3D-reconstruction by least-squares on the non-rectangular reconstruction domain corresponding to the
images of the bust. (f) 3D-reconstruction using the Mumford-Shah approach, on the whole grid. When discontinuities are
handled, it is possible to perform photometric stereo without prior segmentation of the object.

ods can already compete with the compression rate of 4. Aubert, G., Kornprobst, P.: Mathematical Problems in
the well-known JPEG 2000 algorithm [45]. We believe Image Processing, Applied Mathematical Sciences, vol.
147. Springer-Verlag (2002) 2, 8, 13, 18
that the proposed edge-preserving framework may yield
5. Badri, H., Yahia, H., Aboutajdine, D.: Robust Surface
even better results. Reconstruction via Triple Sparsity. In: Proceedings of
Eventually, some of the research directions already the IEEE Conference on Computer Vision and Pattern
mentioned in the conclusion section of our survey pa- Recognition, pp. 2291–2298. Columbus, USA (2014) 14,
per [48] were ignored in this second paper, but they 18
6. Bähr, M., Breuß, M., Quéau, Y., Bouroujerdi, A.S.,
remain of important interest. One of the most appeal-
Durou, J.D.: Fast and accurate surface normal integra-
ing examples is multi-view normal field integration [15]. tion on non-rectangular domains. Computational Visual
Indeed, discontinuities represent a difficulty in our case Media 3, 107–129 (2017) 4, 6
because they are induced by occlusions, yet more infor- 7. Beck, A., Teboulle, M.: A Fast Iterative Shrinkage-
mation would be obtained near the occluding contours Thresholding Algorithm for Linear Inverse Problems.
SIAM Journal on Imaging Sciences 2(1), 183–202 (2009)
by using additional views. 13
8. Belhachmi, B., Bucur, D., Burgeth, B., Weickert, J.: How
to Choose Interpolation Data in Images. SIAM Journal
References on Applied Mathematics 70(1), 333–352 (2009) 21
9. Boyd, S., Parikh, N., Chu, E., Peleato, B., Eckstein, J.:
1. Agrawal, A., Raskar, R., Chellappa, R.: What Is the
Distributed Optimization and Statistical Learning via
Range of Surface Reconstructions from a Gradient Field?
the Alternating Direction Method of Multipliers. Founda-
In: Proceedings of the 9th European Conference on Com-
tions and Trends in Machine Learning 3(1), 1–122 (2011)
puter Vision (volume I), Lecture Notes in Computer Sci-
12
ence, vol. 3951, pp. 578–591. Graz, Austria (2006) 6, 8,
16, 18 10. Bredies, K., Holler, M.: A TGV-Based Framework for
2. Ambrosio, L., Tortorelli, V.M.: Approximation of Func- Variational Image Decompression, Zooming, and Recon-
tionals Depending on Jumps by Elliptic Functionals via struction. Part I: Analytics. SIAM Journal on Imaging
Γ -convergence. Communications in Pure and Applied Sciences 8(4), 2814–2850 (2015) 13, 21
Mathematics 43, 999–1036 (1990) 19 11. Catté, F., Lions, P.L., Morel, J.M., Coll, T.: Image Selec-
3. Attouch, H., Buttazzo, G., Michaille, G.: Variational tive Smoothing and Edge Detection by Nonlinear Diffu-
analysis in Sobolev and BV spaces: applications to PDEs sion. SIAM Journal on Numerical Analysis 29(1), 182–
and optimization. SIAM (2014) 4 193 (1992) 2
Variational Methods for Normal Integration 23

(a) (b) (c)

Fig. 16 Application to image compression/image editing. (a) Reference image. (b) Control points (where the RGB-values and
their gradients are kept). (c) Restored image obtained by considering the proposed Mumford-Shah integrator as a piecewise-
constant interpolation method. A reasonable piecewise constant restoration of the initial image can be obtained from as few
as 10% of the initial information.

12. Chambolle, A.: Image Segmentation by Variational 22. Geman, D., Reynolds, G.: Constrained Restoration and
Methods: Mumford and Shah Functional and the Discrete Recovery of Discontinuities. IEEE Transactions on Pat-
Approximation. SIAM Journal of Applied Mathematics tern Analysis and Machine Intelligence 14(3), 367–383
55(3), 827–863 (1995) 19 (1992) 2
13. Chambolle, A., Pock, T.: A First-Order Primal-Dual Al- 23. Goldstein, T., O’Donoghue, B., Setzer, S., Baraniuk, R.:
gorithm for Convex Problems with Applications to Imag- Fast Alternating Direction Optimization Methods. SIAM
ing. Journal of Mathematical Imaging and Vision 40(1), Journal on Imaging Sciences 7(3), 1588–1623 (2014) 2,
120–145 (2010) 2, 13 13
14. Chan, T.F., Vese, L.A.: Active Contours Without Edges. 24. Goldstein, T., Osher, S.: The split Bregman method for
IEEE Transactions on Image Processing 10(2), 266–277 L1-regularized problems. SIAM Journal on Imaging Sci-
(2001) 19 ences 2(2), 323–343 (2009) 12
15. Chang, J.Y., Lee, K.M., Lee, S.U.: Multiview Normal 25. Haque, S.M., Chatterjee, A., Govindu, V.M.: High Qual-
Field Integration Using Level Set Methods. In: Pro- ity Photometric Reconstruction Using a Depth Camera.
ceedings of the IEEE Conference on Computer Vision In: Proceedings of the IEEE Conference on Computer Vi-
and Pattern Recognition, Workshop on Beyond Multi- sion and Pattern Recognition, pp. 2283–2290. Columbus,
view Geometry: Robust Estimation and Organization USA. (2014) 3
of Shapes from Multiple Cues. Minneapolis, Minnesota,
26. Harker, M., O’Leary, P.: Regularized Reconstruction of
USA (2007) 22
a Surface from its Measured Gradient Field. Journal of
16. Charbonnier, P., Blanc-Féraud, L., Aubert, G., Bar- Mathematical Imaging and Vision 51(1), 46–70 (2015) 4,
laud, M.: Deterministic Edge-Preserving Regularization 6, 8, 9, 10, 18, 21
in Computed Imaging. IEEE Transactions on Image Pro-
27. Hayya, J., Armstrong, D., Gressis, N.: A note on the
cessing 6(2), 298–311 (1997) 2
ratio of two normally distributed variables. Management
17. Chen, S.S., Donoho, D.L., Saunders, M.A.: Atomic De-
Science 21(11), 1338–1341 (1975) 17
composition by Basis Pursuit. SIAM Journal on Scientific
Computing 20(1), 33–61 (1998) 13 28. Hoeltgen, L., Setzer, S., Weickert, J.: An Optimal Control
Approach to Find Sparse Data for Laplace Interpolation.
18. Du, Z., Robles-Kelly, A., Lu, F.: Robust Surface Re-
In: Proceedings of the 9th International Workshop on En-
construction from Gradient Field Using the L1 Norm.
ergy Minimization Methods in Computer Vision and Pat-
In: Proceedings of the 9th Biennial Conference of the
tern Recognition, Lecture Notes in Computer Science,
Australian Pattern Recognition Society on Digital Image
vol. 8081, pp. 151–164. Lund, Sweden (2013) 21
Computing Techniques and Applications, pp. 203–209.
Glenelg, Australia (2007) 12 29. Horn, B.K.P., Brooks, M.J.: The Variational Approach to
19. Durou, J.D., Aujol, J.F., Courteille, F.: Integration of a Shape From Shading. Computer Vision, Graphics, and
Normal Field in the Presence of Discontinuities. In: Pro- Image Processing 33(2), 174–208 (1986) 4
ceedings of the 7th International Workshop on Energy 30. Horovitz, I., Kiryati, N.: Depth from Gradient Fields and
Minimization Methods in Computer Vision and Pattern Control Points: Bias Correction in Photometric Stereo.
Recognition, Lecture Notes in Computer Science, vol. Image and Vision Computing 22(9), 681–694 (2004) 3
5681, pp. 261–273. Bonn, Germany (2009) 1, 2, 14 31. Ikehata, S., Wipf, D., Matsushita, Y., Aizawa, K.: Pho-
20. Durou, J.D., Courteille, F.: Integration of a Normal Field tometric Stereo Using Sparse Bayesian Regression for
without Boundary Condition. In: Proceedings of the 11th General Diffuse Surfaces. IEEE Transactions on Pat-
IEEE International Conference on Computer Vision, 1st tern Analysis and Machine Intelligence 36(9), 1816–1831
Workshop on Photometric Analysis for Computer Vision. (2014) 11
Rio de Janeiro, Brazil (2007) 1, 2, 3, 6 32. Kadambi, A., Taamazyan, V., Shi, B., Raskar, R.: Polar-
21. Gabay, D., Mercier, B.: A dual algorithm for the solu- ized 3D: High-Quality Depth Sensing With Polarization
tion of nonlinear variational problems via finite element Cues. In: Proceedings of the 15th IEEE International
approximation. Computers & Mathematics with Appli- Conference on Computer Vision, pp. 3370–3378. Santi-
cations 2(1), 17 – 40 (1976) 12 ago, Chili (2015) 3
24 Yvain Quéau et al.

33. Kimmel, R., Yavneh, I.: An Algebraic Multigrid Ap- 50. Saracchini, R.F.V., Stolfi, J., Leitão, H.C.G., Atkinson,
proach for Image Analysis. SIAM Journal on Scientific G.A., Smith, M.L.: A Robust Multi-Scale Integration
Computing 24(4), 1218–1231 (2003) 3 Method to Obtain the Depth From Gradient Maps. Com-
34. Kornprobst, P., Aubert, G.: Image Sequence Analysis via puter Vision and Image Understanding 116(8), 882–895
Partial Differential Equations. Journal of Mathematical (2012) 11
Imaging and Vision 11(1), 5–26 (1999) 2 51. Shefi, R., Teboulle, M.: Rate of Convergence Analysis of
35. Koutis, I., Miller, G.L., Peng, R.: A Nearly-m log n Time Decomposition Methods Based on the Proximal Method
Solver for SDD Linear Systems. In: Proceedings of the of Multipliers for Convex Minimization. SIAM Journal
IEEE Annual Symposium on Foundations of Computer on Optimization 24(1), 269–297 (2014) 13
Science, pp. 590–598. Palm Springs, USA (2011) 6 52. Shi, B., Wu, Z., Mo, Z., Duan, D., Yeung, S.K., Tan, P.: A
36. Lanza, A., Morigi, S., Sgallari, F.: Convex Image Denois- Benchmark Dataset and Evaluation for Non-Lambertian
ing via Non-convex Regularization with Parameter Selec- and Uncalibrated Photometric Stereo. In: Proceedings of
tion. Journal of Mathematical Imaging and Vision 56(2), the IEEE Conference on Computer Vision and Pattern
195–220 (2016) 2, 14 Recognition. Las Vegas, USA (2016) 17, 18
37. Mumford, D., Shah, J.: Optimal Approximations by 53. Simchony, T., Chellappa, R., Shao, M.: Direct Analytical
Piecewise Smooth Functions and Associated Variational Methods for Solving Poisson Equations in Computer Vi-
Problems. Communications in Pure and Applied Math- sion Problems. IEEE Transactions on Pattern Analysis
ematics 42(5), 577–685 (1989) 2, 11, 18 and Machine Intelligence 12(5), 435–446 (1990) 6, 8, 9,
38. Nikolova, M.: Local Strong Homogeneity of a Regular- 10, 18
ized Estimator. SIAM Journal on Applied Mathematics 54. Weickert, J.: Anisotropic diffusion in image processing.
61(2), 633–658 (2000) 2, 13 Teubner Stuttgart (1998) 16
39. Noakes, L., Kozera, R.: Nonlinearities and Noise Reduc- 55. Woodham, R.J.: Photometric Method for Determining
tion in 3-Source Photometric Stereo. Journal of Mathe- Surface Orientation from Multiple Images. Optical En-
matical Imaging and Vision 18(2), 119–127 (2003) 17 gineering 19(1), 139–144 (1980) 4, 16, 17, 18, 21
40. Ochs, P., Brox, T., Pock, T.: iPiasco: Inertial proximal
algorithm for strongly convex optimization. Journal of
Mathematical Imaging and Vision 53(2), 171–181 (2015)
2
41. Ochs, P., Chen, Y., Brox, T., Pock, T.: iPiano: Inertial
Proximal Algorithm for Nonconvex Optimization. SIAM
Journal on Imaging Sciences 7(2), 1388–1419 (2014) 2,
14
42. Or-el, R., Rosman, G., Wetzler, A., Kimmel, R., Bruck-
stein, A.M.: RGBD-Fusion: Real-Time High Precision
Depth Recovery. In: Proceedings of the IEEE Confer-
ence on Computer Vision and Pattern Recognition, pp.
5407–5416. Boston, USA. (2015) 3
43. Pérez, P., Gangnet, M., Blake, A.: Poisson image editing.
ACM Transactions on Graphics 22(3), 313–318 (2003) 21
44. Perona, P., Malik, J.: Scale-space and Edge Detection us-
ing Anisotropic Diffusion. IEEE Transactions on Pattern
Analysis and Machine Intelligence 12(7), 629–639 (1990)
2, 11, 16, 17
45. Peter, P., Hoffmann, S., Nedwed, F., Hoeltgen, L., We-
ickert, J.: Evaluating the true potential of diffusion-based
inpainting in a compression context. Signal Processing:
Image Communication 46, 40 – 53 (2016) 21, 22
46. Pock, T., Cremers, D., Bischof, H., Chambolle, A.: An
algorithm for minimizing the Mumford-Shah functional.
In: Proceedings of the 12th IEEE International Confer-
ence on Computer Vision, pp. 1133–1140. Kyoto, Japan
(2009) 19
47. Quéau, Y., Durou, J.D.: Edge-Preserving Integration of
a Normal Field: Weighted Least Squares, TV and L1 Ap-
proaches. In: Proceedings of the 5th International Con-
ference on Scale Space and Variational Methods in Com-
puter Vision, Lecture Notes in Computer Science, vol.
9087, pp. 576–588. Lège Cap-Ferret, France (2015) 1, 2,
11, 12, 16
48. Quéau, Y., Durou, J.D., Aujol, J.F.: Normal Integration:
A Survey. Journal of Mathematical Imaging and Vision
(2017). (submitted, preprint available at [Link]
[Link]/hal-01334349/) 1, 3, 10, 21, 22
49. Rudin, L.I., Osher, S., Fatemi, E.: Nonlinear total varia-
tion based noise removal algorithms. Physica D: Nonlin-
ear Phenomena 60(1-4), 259–268 (1992) 12