A unified approach to the well-posedness of some non-Lambertian models in Shape-from-Shading theory Fabio Camilli∗ Silvia Tozza† arXiv:1603.06336v3 [math.AP] 4 Oct 2016 October 11, 2018 Abstract In this paper we show that the introduction of an attenuation factor in the brightness equations relative to various perspective Shape-from-Shading models allows to make the corresponding differential problems well-posed. We propose a unified approach based on the theory of viscosity solutions and we show that the brightness equations with the attenuation term admit a unique viscosity solution. We also discuss in detail the possible boundary conditions that we can use for the Hamilton-Jacobi equations associated to these models. Keywords: Shape-from-Shading, perspective model, viscosity solution, maximum principle, non-Lambertian models, stationary Hamilton-Jacobi equations. 1 Introduction Shape-from-Shading (SfS) is the problem to compute a three dimensional shape of a surface from a single gray-value image of it. The SfS model was formally introduced by Horn who first formulated the problem as a nonlinear first order partial differential equation (PDE) of Hamilton-Jacobi (HJ) type. In general, this problem is described by the following so-called image irradiance equation introduced by Bruss [8] I(x) = R(N (x)). (1) In this equation, the normalized brightness I(x) of the given gray-value image is put in relation with the reflectance map R(N (x)) giving the value of the light reflection on the surface as a function of its orientation (i.e., of the normal N (x)). Depending on how we describe the function R, different reflection models are determined. The most studied model is the Lambertian one, a view independent model that depends only on the incident angle between the normal and the light source direction ω. Considering this model under an orthographic projection, that is when a single light source is placed at infinity in the direction of the unit vector ω and assuming uniform albedo equal to 1 (i.e. the light is completely reflected by the surface) the irradiance equation (1) becomes p (2) I(x) 1 + |∇u(x)|2 + ω · (∇u(x), −1) = 0, where u is the surface height (our unknown). In general, nonlinear PDEs such as (2) do not admit smooth solutions and the proper notion of weak solution is the viscosity one (see [4, 9, 3]). Viscosity solutions in the framework of SfS were considered for the first time in [24] and since then many papers have been devoted to develop a proper existence and uniqueness theory for (2) and related models. It is worthwhile to note that the orthographic Lambertian model suffers of two main drawbacks: ∗ Dip. di Scienze di Base e Applicate per l’Ingegneria, “Sapienza” Università di Roma, via Scarpa 16, 00161 Rome, Italy, (e-mail:camilli@sbai.uniroma1.it) † Dipartimento di Matematica, “Sapienza” Università di Roma, P.le Aldo Moro, 5 - 00185 Rome, Italy (e-mail: tozza@mat.uniroma1.it) 1 (i) From a modeling point of view, the Lambertian model does not take into account the viewer direction and supposes to deal with smooth objects, assumptions not suitable to deal with realworld data. Moreover, except for few application fields such as the astronomical one, the light source and the camera are not far away from the surface (e.g. applications in medicine, security, etc.). (ii) The image irradiance equation (2), even if completed with boundary conditions, does not admit in general a unique solution. In fact, due to the well-known concave/convex ambiguity in presence of singular points, i.e. points where the light direction is parallel to the normal to the surface, there exists an infinite family of viscosity solutions to equation (2) which are all in between a minimal and a maximal solution. Concerning the first issue, in order to increase the complexity of the model, various extensions have been proposed. In [22, 17], the orthographic Lambertian model is replaced by a perspective Lambertian one (pinhole camera and light source at the optical center). In [16], a more general setup, which combines a spherical surface parametrization with the non-Lambertian Oren-Nayar reflectance model, is considered obtaining a robust approach that allows to deal with an arbitrary position of the light source. Moreover, models for non-Lambertian surfaces (like e.g. [21, 6, 20]) have been proposed in [1, 2, 26, 16]. All these models, when formulated in terms of a differential equation, do not resolve the concave/convex ambiguity (see [25] for an analysis on them). In order to solve the non-uniqueness issue in presence of singular points, in [23, 22] the authors introduce in the perspective Lambertian model an attenuation factor which takes into account the distance between the light source and the surface. This new term makes the model more reliable and thanks to it the notion of singular points does not make sense anymore (even if some ambiguities in the model still appear [7]). They show that a Hamilton-Jacobi equation, obtained via an appropriate change of variable in the image irradiance equation, admits a unique viscosity solution and, therefore, the uniqueness of the solution for the original equation follows. In this paper, starting from the idea in [22] related to the Lambertian model, we consider non-Lambertian models (Oren-Nayar [20], Phong [21], Blinn-Phong [6]) and we show that the introduction of an attenuation factor in the perspective case allows to overcome the concave/convex ambiguity for all the corresponding brightness equations and to obtain therefore the existence of at most one viscosity solution. The uniqueness property is obtained by checking the assumptions of the Maximum Principle for discontinuous viscosity solutions (see [4]). In this framework, boundary conditions play an important role, hence we discuss the various type of them which is possible to prescribe for the different models. Interesting enough, all the models (included the Lambertian one) can be studied in a unified approach which allows to check the assumptions of the Maximum Principle in a rather straightforward way. Moreover, this approach can also provide a setup to study other SfS perspective models (see [2]). Many papers have been devoted to the approximation of the various irradiance equations arising in SfS theory (see the surveys [27, 12]), but only very few of them discuss the convergence of the corresponding methods and most of the results are of experimental type. This fact is essentially due to the lack of a proper characterization of the viscosity solution of the limit equation. The Maximum Principle for discontinuous viscosity solutions we prove here is a key ingredient to show the convergence of numerical schemes for SfS equations via the by now classical method proposed in [5] (we will discuss this point elsewhere). The paper is organized as follows: in Section 2, we briefly describe the reflectance models and the derivation of the PDEs associated to them, arriving to a unified formulation with common properties shared by the different models. In Section 3, we focus our attention on viscosity solutions and the role of the various possible boundary conditions for HJ equations. After that, we move to prove the well-posedness of the Perspective SfS problem for the non-Lambertian models in Section 4. Finally, we report final comments in Section 5 and in Appendix A the proofs of two lemmas used for the unified analysis. 2 Reflectance Models and PDE formulation In this section we will briefly recall the reflectance models we will consider and we will describe the derivation of the corresponding HJ equations, arriving to a unified formulation of them. For all the models, we consider as setup the perspective camera projection with one light source 2 located at the optical center of the camera and an attenuation term of the illumination due to the distance between the light source and the surface. Let Ω be an open subset of R2 representing the image domain. Since the CCD sensors have finite size, let us assume that Ω is bounded. As in [23, 22, 1, 26], we represent the scene by a surface parametrized by the function S : Ω → R3 defined by S(x) = p f u(x) |x|2 + f 2 (x, −f ) where f > 0 is the focal length of the camera and u(x) is the unknown height of the surface. The normal vector at the point x := (x1 , x2 ) of S is given by
f u(x) f u(x) x, ∇u(x) · x + 2 f 2 2 2 |x| + f |x| + f n(x) = Sx1 × Sx2 = f ∇u(x) − (3) and the unit vector in the direction of the illumination is (−x, f ) ω(S(x)) = p . |x|2 + f 2 The unit normal vector N (x) will be N (x) = n(x) |n(x)| , where s |n(x)| = 2.1 f 2 |∇u(x)|2 + (∇u(x) · x)2 + u(x)2 f2 . (f 2 + |x|2 ) (4) Lambertian model (L–model) The Lambertian model is the most common reflectance model ([22, 17, 14, 15]). The brightness equation for this model is described by I(x) = γD (x) cos(θi ) r2 (5) where I(x) denotes the image intensity, θi represents the angle between the unit normal to the surface N (x) and the light source direction ω(S(x)), γD (x) is the diffuse albedo which describes the physical properties of the surface reflection, and 1/r2 is the attenuation term of the illumination due to the distance between the light source and the surface [1]. Hence, for a Lambertian surface, the measured light in the image depends only on the scalar product between N and ω, independently on the position of the viewer. In the sequel, we suppose uniform albedo and we put γD (x) = 1, that is all the points of the surface reflect completely the light that hits them (we will not consider absorption phenomena). So in this case the surface of the object has uniform properties of light reflection (see Fig. 1). Under a perspective camera projection with a single light source located at the optical center of the camera and considering the attenuation term 1/r2 , with r = f u(x), we can associate to (5) the following HJ equation I(x) = 1 f 2 u(x)2 u(x)Q(x) p 2 2 f |∇u(x)| + (∇u(x) · x)2 + u(x)2 Q(x)2 where Q(x) := p f f2 + |x|2 . (6) (7) By defining W (x, p) := f 2 |p|2 + (p · x)2 , Q2 (8) we can rewrite Eq. (6) as I(x) = 1 f 2 u(x)2 u(x) p W (x, ∇u(x)) + u(x)2 3 . Figure 1: A Lambertian surface diffuses the incident radiation independently from the angle between the surface normal N and the direction of the observer. Two different observers V1 and V2 do not detect any difference in radiance. The radiance is a function only of the angle θi between ω and N . Assuming that the surface S is visible, i.e. in front of the optical center, then u is a strictly positive function. Hence, by the change of variable v(x) = ln u(x), we get the following HJ equation in the new variable v (as in [23, 22]) p − e−2v(x) + I(x)f 2 W (x, ∇v(x)) + 1 = 0, (9) to which we associate the Hamiltonian H L (x, r, p) = −e−2r + I(x)f 2 2.2 p W (x, p) + 1 . (10) Oren-Nayar model (ON–model) The Oren-Nayar reflectance model [20, 19] represents a rough surface as an aggregation of Vshaped cavities, each with Lambertian reflectance properties. Assuming that there is a linear relation between the irradiance of the image and the image intensity, the brightness equation for the ON–model is given by I(x) = cos(θi )(A + B sin(α) tan(β)M (ϕi , ϕr )). (11) In this model, θi represents the angle between the unit normal to the surface N and the light source direction ω, θr stands for the angle between N and the observer direction V , ϕi is the angle between the projection of the light source direction ω and the x1 axis onto the (x1 , x2 )plane, ϕr denotes the angle between the projection of the observer direction V and the x1 axis onto the (x1 , x2 )-plane, the two variables α and β are given by α = max {θi , θr } and β = min {θi , θr } . The nonnegative constants A and B depend on the statistics of the cavities via the roughness parameter σ. We set σ ∈ [0, π/2), interpreting σ as the slope of the cavities, and A = 1 − 0.5 σ 2 (σ 2 + 0.33)−1 B = 0.45σ 2 (σ 2 + 0.09)−1 M (ϕi , ϕr ) = max{0, cos(ϕr − ϕi )}. Let us consider a perspective camera projection, with only one light source located at the optical center of the camera. As a consequence of this setup, we define θ := θi = θr = α = β. Adding the attenuation term 1/r2 of the illumination due to the distance between the light source and the surface as proposed in [1], the brightness equation becomes I(x) = A cos(θ) + B sin2 (θ) . r2 4 (12) Since cos(θ) = N (x) · ω(S(x)), hence sin2 (θ) = 1 − (N (x) · ω(S(x)))2 , and r = f u(x), the brightness equation (12) becomes Au(x)Q(x) p 2 2 f |∇u| + (∇u · x)2 + u(x)2 Q(x)2   u(x)2 Q(x)2 +B 1 − 2 f |∇u|2 + (∇u · x)2 + u(x)2 Q(x)2 I(x)f 2 u(x)2 = (13) with Q(x) defined as in (7). By defining W as in (8), we can rewrite Eq. (13) as I(x)f 2 u(x)2 − p Au(x) W (x, ∇u(x)) + u(x)2 − B(1 − u(x)2 )=0 W (x, ∇u(x)) + u(x)2 or, equivalently, p I(x)f 2 u(x)2 (W (x, ∇u(x)) + u(x)2 ) − Au(x)( W (x, ∇u(x)) + u(x)2 ) − BW (x, ∇u(x)) = 0. Since u(x) > 0, ∀x ∈ Ω, by the change of variable v(x) = ln u(x) we arrive to the following PDE in v p If 2 e2v(x) (W (x, ∇v(x)) + 1) − A W (x, ∇v(x)) + 1 − BW (x, ∇v(x)) = 0, which can be written also as follow (Cf. Eq. (10) in [1]): (W (x, ∇v(x)) + 1) = 0, − e−2v(x) + I(x)f 2 p A W (x, ∇v(x)) + 1 + BW (x, ∇v(x)) to which we associate the following Hamiltonian " H ON (x, r, p) := −e −2r (W (x, p) + 1) 2 + I(x) f A p W (x, p) + 1 + BW (x, p) (14) # . (15) Note that the quantity W is the same one defined in the L–model, see (8), and this will be useful for the unified approach to the different models. 2.3 Phong model (PH–model) The PH–model is an empirical model that was developed by Phong [21]. This model introduces a specular component to the brightness function I(x) which can be described in general as kA IA (x) + kD ID (x) + kS IS (x) where IA (x), ID (x) and IS (x) are the ambient, diffuse and specular light components, respectively, and kA , kD and kS indicate the percentages of ambient, diffuse and specular components, respectively, such that their sum is equal to 1. The specular light component IS (x) is given as a power of the cosine of the angle θs between the unit vectors V and R(x) with R representing the reflection of the light ω on the surface. Then, the brightness equation for the PH–model is I(x) = kA IA (x) + kD γD (x)(cos θi ) + kS γS (x)(cos θs )α , (16) where α expresses the specular reflection characteristics of a material, γD (x) and γS (x) represent the diffuse and specular albedo, respectively (we will set them to 1 in the following). Using the same setup adopted for the previous models, that is a perspective camera projection with only one light source located at the optical center of the camera, we obtain that the view direction and the light source direction are the same, hence θs = 2θi (see [26]). By adding the light attenuation term as done before, the brightness equation (16) becomes I(x) = kA IA (x) + 1 kD (N (x) · ω(S(x))) + kS (2(N (x) · ω(S(x)))2 − 1)α , r2 (17) since cos θi = N (x) · ω(S(x)) and cos θs = cos 2θi = 2(cos θi )2 − 1 = 2(N (x) · ω(S(x)))2 − 1. 5 Following [26], we derive a HJ equation associated to (17). If we use |n(x)| as defined in (4) and Q(x) as in (7), then (17) is equivalent to I(x) = kA IA (x) +   u(x)Q(x) 1 u2 Q(x)2 α . k − 1) + k (2 D S f 2 u(x)2 |n(x)| |n(x)|2 (18) Since u(x) is a strictly positive function, by the change of variable v(x) = ln u(x) we arrive to the following HJ equation in v (I(x) − kA IA (x))f 2 V BW W (x, ∇v(x)) − kD e−2v(x) Q α kS W V BW (x, ∇v(x)) −2v(x)  2Q(x)2 − 1 = 0, − e Q(x) W V BW (x, ∇v(x))2 (19) where W V BW (x, p) := p f 2 |p|2 + (p · x)2 + Q(x)2 , as defined in [26]. To our purpose, it is useful to write (19) with the same structure and functions that appear in the L-model and in the ON-model. Hence, we rewrite (19) as follow p W (x, ∇v(x)) + 1 −2v(x) 2  α = 0, −e + (I(x) − kA IA (x))f (20) p (x,∇v(x)) kD + kS W (x, ∇v(x)) + 1 1−W 1+W (x,∇v(x)) where W (x, ∇v(x)) is defined as in (8). We associate to (20) the following Hamiltonian p W (x, p) + 1 PH −2r 2  α . H (x, r, p) := −e + (I(x) − kA IA (x))f p (x,p) kD + kS W (x, p) + 1 1−W 1+W (x,p) (21) From now on we assume that α is an integer and α ≥ 1 (in [21], α ∈ [1, 10]). Note that H P H (x, r, p) is not well  defined if α is irrational. On the other side, as stated also in Vogel et al. 1−W (x,p) 1+W (x,p) represents the cosine of the specular term, in the Phong model   (x,p) this cosine is replaced by zero if cos θs = 1−W 1+W (x,p) < 0 and we can rewrite (21) as follow: [26], since the term H P H (x, r, p) := p  W (x, p) + 1  −2r 2  + (I(x) − kA IA (x))f   −e kD p W (x, p) + 1   α  −e−2r + (I(x) − kA IA (x))f 2 p  (x,p)  kD + kS W (x, p) + 1 1−W 1+W (x,p) if W (x, p) ≥ 1, (22) if 0 ≤ W (x, p) < 1. Remark 2.1. Note that if kS = 0, i.e. if there is no contribution of the specular part, then   p I(x) − kA IA (x) H P H (x, r, p) = −e−2r + f 2 W (x, p) + 1. kD Hence, in this case the PH–model reduces to the L–model, up to a constant. 2.4 Blinn-Phong model (BP–model) A modification of the PH–model has been proposed by Blinn [6] introducing an intermediate vector H, which bisects the angle between the unit vectors ω and V . If the surface is a perfect mirror, the light reaches the eye only if the surface normal N is pointed halfway between the light source direction ω and the eye direction V . We will denote the direction of maximum ω+V highlight by H = |ω+V | . For less than perfect mirrors, the specular component falls off slowly as the normal direction moves away from the specular direction. The cosine of the angle between N and H is used as a measure of the distance of a particular surface from the maximum specular direction. The degree of sharpness of the highlights is adjusted by taking this cosine to some 6 power (in [6] it is stated that this power is typically 50 or 60). For this model, the brightness equation is I(x) = kA IA (x) + kD γD (x)(cos θi ) + kS γS (x)(cos δ)c , (23) where c expresses a measure of shininess of the surface (we suppose c ≥ 1) and δ is the angle between the unit normal N and H. As for the PH–model, γD (x) and γS (x) represent the diffuse and specular albedo, respectively, and in what follows we will set them to 1. Using the same setup adopted for the previous models, we obtain that the view direction and the light source direction are the same. By adding the light attenuation factor as done before, the brightness equation (23) becomes I(x) = kA IA (x) + 1 kD (N (x) · ω(S(x))) + kS (N (x) · ω(S(x)))c , r2 since cos θi = N (x) · ω(S(x)) and cos δ can be expressed as N (x) · ω(S(x)) under this setup. By using (3) and (4) for the definition of the unit normal N (x), we arrive to the following HJ equation  u(x)Q(x)  u(x)Q(x) c  1 + kS I(x) = kA IA (x) + 2 kD , 2 f u(x) |n(x)| |n(x)| with |n(x)| defined as in (4) and Q(x) as in (7). By using the change of variable v(x) = ln u(x), as done for the previous models, we obtain the following HJ equation in v (I(x) − kA IA (x))f 2 e2v = kD p 1 W (x, ∇v(x)) + 1  + kS p c 1 W (x, ∇v(x)) + 1 (24) with the function W (x, ∇v(x)) defined as in (8). For our purpose, we will rewrite this equation as (W (x, ∇v(x)) + 1)c/2 − e−2v(x) + (I(x) − kA IA (x))f 2 = 0, (25) c−1 kD (W (x, ∇v(x)) + 1) 2 + kS to which we associate the Hamiltonian H BP (x, r, p) := −e−2r + (I(x) − kA IA (x))f 2 (W (x, p) + 1)c/2 kD (W (x, p) + 1) c−1 2 + kS . (26) Remark 2.2. For the Blinn-Phong model holds an observation similar to the one reported for the PH–model in Remark 2.1 , i.e. for kS = 0 it reduces to the L-model, up to a constant. 2.5 Unified formulation and analysis As we have seen before, the four different reflectance models considered in a common setup lead to the study of a HJ equation of the form H(x, v(x), ∇v(x)) = 0, ∀x ∈ Ω, (27) where the Hamiltonian H is given respectively by (10), (15), (22), (26) and v(x) = ln u(x), with u(x) the unknown height of the surface. We highlight some properties of the general Hamiltonian H which are common for all the four cases. First of all, note that (27) can be rewritten as − e−2v + HM (x, ∇v(x)) = 0 (28) where M is the acronym of the models (M = L, ON, P H, BP ). For W (x, p) defined as in (8), HM (x, p) corresponds for each model to the following: Lambertian model HL (x, p) := I(x)f 2 FL (W (x, p)), √ FL (r) = r + 1. 7 (29) (30) Oren-Nayar model HON (x, p) = I(x)f 2 FON (W (x, p)), r+1 . FON (r) = √ A r + 1 + Br (31) (32) Phong model HP H (x, p) = (I(x) − kA IA (x))f 2 FP H (W (x, p)),  √ r+1    kD FP H (r) = (r + 1)α+1/2    kD (r + 1)α + kS (r + 1)1/2 (1 − r)α (33) if r ≥ 1 (34) if 0 ≤ r < 1. Blinn-Phong model HBP (x, p) = (I(x) − kA IA (x))f 2 FBP (W (x, p)), FBP (r) = (r + 1) kD (r + 1) (35) c/2 c−1 2 + kS . (36) We give some bounds on W and the various functions FM above defined which will be useful for our desired unified analysis. The proofs of these results are reported at the end of the paper in Appendix A. Lemma 2.1. The following bounds hold: f 2 |p|2 ≤ W (x, p) ≤ C|p|2 |Dx W (x, p)| ≤ C|p| 2 (37) (38) |Dp W (x, p)| ≤ C|p| (39) with C independent of (x, p) ∈ Ω × R2 . Lemma 2.2. (i) For F = FL , FP H , FBP , the function F : [0, +∞) → R is smooth, positive, strictly increasing and satisfies lim F (r) = +∞, (40) C . r+1 (41) r→+∞ 0 < F 0 (r) ≤ √ (ii) The function FON : [0, +∞) → R is smooth, positive, strictly increasing if A/2 > B and satisfies 1 , B (42) C . r+1 (43) lim FON (r) = r→+∞ 0 0 < FON (r) ≤ √ 3 Viscosity solutions of Hamilton-Jacobi equations In view of the study of the general Hamiltonian (27) for the various reflectance models introduced in Section 2, we recall some basic properties of the theory of viscosity solutions emphasizing the role of the boundary conditions (we refer to [4] for a nice introduction to this theory). In general, a Hamilton-Jacobi equation such as (27) does not admit a classical solution and the right notion of weak solution is the one of viscosity solution. Viscosity solutions are typically Lipschitz continuous and they can be obtained as the limit of regular solutions to second order 8 problems via the so-called vanishing viscosity method. Moreover, if a classical solution to the HJ equation exists, it coincides with the viscosity solution. We introduce the definition of viscosity subsolution, supersolution and solution. It is useful to have the following notations for a generic set O ⊂ R2 USC(O) = {upper semicontinuous functions u : O → R} LSC(O) = {lower semicontinuous functions u : O → R} Definition 3.1. A viscosity subsolution of (27) on Ø is a function u ∈ USC(Ø) such that for any φ ∈ C 1 (Ø), if x0 ∈ Ø is a local maximum point of u − φ, one has H(x0 , u(x0 ), ∇φ(x0 )) ≤ 0. A viscosity supersolution of (27) on Ø is a function u ∈ LSC(Ø) such that for any φ ∈ C 1 (Ø), if x0 ∈ Ø is a local minimum point of u − φ, one has H(x0 , u(x0 ), ∇φ(x0 )) ≥ 0. Finally, u is a viscosity solution of (27) if it is both a viscosity subsolution and a viscosity supersolution. Note that a viscosity solution of (27) is a continuous function in Ø. 3.1 Strong boundary conditions In this framework we assume to know the value of the function u on the boundary of the image domain Ø, hence we consider the Dirichlet problem:  H(x, u(x), ∇u(x)) = 0, ∀x ∈ Ω, (44) u(x) = g(x), ∀x ∈ ∂Ω, where g is a real continuous function defined on ∂Ω and the boundary condition is considered in pointwise sense (boundary condition in strong sense). The following theorem ensures the uniqueness of the continuous viscosity solution to (44) (see [9, Theorem 3.3] for more details). Theorem 3.2. Let Ω be a bounded subset of R2 and H : Ω × R × R2 → R be continuous and satisfies for any R > 0, ∃γR > 0 s.t. H(x, u, p) − H(x, v, p) ≥ γR (u − v) for all x ∈ Ω, −R ≤ v ≤ u ≤ R and p ∈ R2 , |H(x, u, p) − H(y, u, p)| ≤ mR (|x − y|(1 + |p|)) for all x, y ∈ Ω, −R ≤ u ≤ R, p ∈ R2 where mR (t) → 0 when t → 0. (H1) (H2) Let u ∈ USC(Ø) (respectively v ∈ LSC(Ø)) be a subsolution (respectively, supersolution) of (27) in Ø and u ≤ v on ∂Ø. Then u ≤ v in Ø. The previous result implies that there exists at most one viscosity solution u ∈ C(Ø) to the Dirichlet problem (44) and we will prove that all the Hamiltonians corresponding to the reflectance models introduced in Section 2 satisfy the assumptions (H1)-(H2). Note that the crucial point for the uniqueness is the assumption (H1), which is in general not satisfied by the SfS models without attenuation factor and/or under a different setup, e.g. considering an orthographic projection instead of a perspective one. However, Theorem 3.2 is not completely satisfying for the following reasons: (i) A Dirichlet boundary condition requires the exact value of the solution at the boundary of the image domain, an additional information which is not always available in the SfS framework, especially in real contexts. (ii) It is well known that boundary conditions for first order PDEs can be typically prescribed only on some part of the boundary (active boundary). Hence, imposing a Dirichlet condition on all the boundary as in (44) could lead to an ill-posed problems, i.e. problems for which there exists no viscosity solution (see [4, cap. IV]) and in this case the corresponding numerical schemes return inconsistent or inexact results. 9 (iii) Many algorithms have been suggested for the SfS problem and a large part of them are based on the approximation of the partial differential equations arising in the various models (see the surveys [27, 12]). Convergence and stability analysis are important points for the reliability of a numerical approximation, but only a very few papers analyze the theoretical properties of various proposed algorithms and most of the results are of empirical nature. It is worthwhile recall that there is a huge literature about the approximation of viscosity solutions (see e.g. the book by Falcone and Ferretti [13]) and a well established technique to prove the convergence of a numerical scheme is the Barles-Souganidis’ approach [5]: besides some natural properties of the scheme (stability, consistency, monotonicity), a key ingredient for this technique is the Maximum Principle for discontinuous viscosity solution, which in particular requires the formulation in a weak (viscosity) sense of the boundary conditions for the PDE to be approximated. Since the Dirichlet boundary condition in (44) is formulated in strong sense, i.e. pointwise, Theorem 3.2 is not sufficient for applying the Barles-Souganidis’ approach and we need a different comparison result. Motivated by the previous remarks, we are going to introduce the notion of boundary conditions in weak (viscosity) sense. 3.2 Weak boundary conditions In this section we consider a more general boundary value problem of the form  H(x, u(x), ∇u(x)) = 0, ∀x ∈ Ω, B(x, u(x), ∇u(x)) = 0, ∀x ∈ ∂Ω. (45) The operator B : ∂Ø × R × R2 → R represents the boundary condition. For example, the Dirichlet boundary condition in (44) arises from the choice B(x, r, p) = r − g(x); the Neumann boundary condition arises if B(x, r, p) = N (x) · p − g(x) where N (x) is the outward unit normal to x ∈ ∂Ø; the so-called state constraints condition corresponds to B(x, r, p) = r − g(x) with g(x) ≡ −∞. We introduce the notion of (discontinuous) viscosity solution for (45). The main difference with Definition 3.1 is that now the boundary condition is incorporated in the definition (boundary condition in weak sense). Definition 3.3. A viscosity subsolution of (45) is a function u ∈ USC(Ø) such that for any φ ∈ C 1 (Ø), if x0 is a local maximum point of u − φ, one has H(x0 , u(x0 ), ∇φ(x0 )) ≤ 0, if x ∈ Ø, min[H(x0 , u(x0 ), ∇φ(x0 )), B(x0 u(x0 ), ∇φ(x0 ))] ≤ 0, if x ∈ ∂Ø. A viscosity supersolution of (45) is a function u ∈ LSC(Ø) such that for any φ ∈ C 1 (Ø), if x0 is a local minimum point of u − φ, one has H(x0 , u(x0 ), ∇φ(x0 )) ≥ 0, if x ∈ Ø, max[H(x0 , u(x0 ), ∇φ(x0 )), B(x0 u(x0 ), ∇φ(x0 ))] ≥ 0, if x ∈ ∂Ø. Finally, u is a viscosity solution of (45) if it is both a viscosity subsolution and a viscosity supersolution. Note that a viscosity solution of (45) is a continuous function in Ø. We can now state the Maximum Principle for discontinuous viscosity solutions. We first consider the case of a Dirichlet boundary condition (see [4, Theorem 4.4]). Theorem 3.4. Let Ω be a bounded subset of R2 with ∂Ø of class W 2,∞ , g : ∂Ø → R a continuous function and H : Ω × R × R2 → R satisfies (H1)-(H2). Moreover, there exists a neighborhood Γ of ∂Ω in R2 such that for all 0 < R < +∞, there exists mR (t) → 0 when t → 0 such that |H(x, u, p) − H(x, u, q)| ≤ mR (|p − q|)), for all x ∈ Γ, −R ≤ u ≤ R, p, q ∈ R2 , 10 (H3) for all 0 < R < +∞, there exists CR > 0 such that H(x, u, p + λn(x)) ≤ 0 ⇒ λ ≤ CR (1 + |p|) for all (x, u, p) ∈ Γ × [−R, R] × R2 , for all R1 , R2 ∈ R+ , H(x, u, p − λn(x)) → +∞ λ → +∞, uniformly for (x, u, p) ∈ Γ × [−R1 , R2 ] × B(0, R2 ). (H4) (H5) Let u ∈ USC(Ø) (respectively v ∈ LSC(Ø)) be a subsolution (respectively, supersolution) of (45) with B(x, r, p) = r − g(x). Then u ≤ v in Ø. Even if the comparison between subsolution and supersolution holds only in Ø and not in all Ø, the previous result implies the uniqueness of the viscosity solution in Ø, i.e. if u1 , u2 ∈ C(Ø) are two viscosity solutions of (45) then u1 = u2 in Ø. Moreover, Theorem 3.4 is the type of result needed to prove the convergence of approximation schemes via the Barles-Souganidis’ method. Concerning the other boundary conditions which are typically considered in the SfS framework, for the validity of the corresponding Maximum Principle it is sufficient to assume (i) for the State Constraints condition B(x, r, p) = r − g(x) with g(x) ≡ −∞, that H is continuous and satisfies (H1)-(H4). Note that in this case the boundary condition reduces to H(x, u, ∇u) ≥ 0 on ∂Ø (46) in viscosity sense and, therefore, it requires no additional information on the solution at the boundary. (ii) for the Neumann condition B(x, r, p) = N (x)·p−g(x), that the function g is continuous on ∂Ø and H is continuous and satisfies (H1)-(H3). In this case the Maximum Principle holds in Ø, i.e. u ≤ v in Ø in Theorem 3.4. 4 Well-posedness of the PSfS problem Aim of this section is to show that the SfS Hamiltonians H L , H ON , H P H , H BP , introduced in Section 2, completed with appropriate (either strong or weak) boundary conditions, give raise to well-posed problems (either (44) or (45)), i.e. problems for which the uniqueness of the viscosity solution holds. In this section we assume that I(x) is Lipschitz continuous in Ω and there exists a neighborhood Γ of ∂Ω such that I(x) ≥ δ > 0 in Γ, where ( I(x) = 4.1 I(x) if H = H L , H ON , I(x) − kA IA (x) if H = H P H , H BP . (47) (48) Validation of assumptions for strong boundary conditions We start now to prove that the assumptions (H1)-(H2) previously defined in Section 3.1 hold for our Hamiltonians. Proposition 4.1. All the SfS Hamiltonians H L , H ON , H P H , H BP defined in (10), (15), (22), (26) satisfy the assumptions (H1)-(H2). Proof. Validation of hypothesis (H1) Recalling that all the Hamiltonians H L , H ON , H P H and H BP can be written in the form H(x, u, p) = −e−2u + HM (x, p) for an appropriate function HM (x, p), see (28), then for u, v such that −R ≤ v ≤ u ≤ R, we have H(x, u, p) − H(x, v, p) = −e−2u + HM (x, p) + e−2v − HM (x, p) = −e−2u + e−2v ≥ 2(u − v)e−2R = γR (u − v). 11 = for some c ∈ [v, u] 2(u − v)e−2c Hence, all the Hamiltonians satisfy assumption (H1) with γR = 2e−2R > 0. Validation of hypothesis (H2) For H(x, u, p) equal to either H L ,H ON , H P H or H BP , by (28) we can write |H(x, u, p) − H(y, u, p)| = |HM (x, p) − HM (y, p)|, (49) where M is the acronym of the models (M = L, ON, P H, BP ). Recalling (47), with I defined as in (48), we can write |HM (x, p) − HM (y, p)| = I(x)f 2 FM (W (x, p)) − I(y)f 2 FM (W (y, p)) ≤ |I(x) − I(y)|f 2 FM (W (x, p)) + f 2 |I(y)||FM (W (x, p)) − FM (W (y, p))| ≤ f 2 kDIk∞ |x − y|FM (W (x, p)) + f 2 kIk∞ |FM (W (x, p)) − FM (W (y, p))|. 0 Since Dx FM (W (x, p)) = FM (W (x, p))Dx W (x, p), we use (37), (38) and (41) for HL , (38), (43) ON PH for H , (38), (41) for H and HBP to obtain (H2) with C depending on kIk∞ and kDIk∞ . As an immediate corollary of the previous Proposition 4.1 and Theorem 3.2, we get Corollary 4.2. The Dirichlet problem (44) for H given by H L , H ON , H P H or H BP admits at most one viscosity solution. 4.2 Validation of assumptions for weak boundary conditions Let us see now the proof of the assumptions (H3)-(H5) introduced in Section 3.2. Proposition 4.3. We have (i) The SfS Hamiltonians H L , H P H , H BP defined in (10), (22), (26), satisfy the assumptions (H3)-(H5). (ii) The SfS Hamiltonian H ON defined in (15) satisfies (H3). Proof. Validation of hypothesis (H3) For H(x, u, p) equal to either H L , H ON , H P H or H BP , we have |Dp H(x, u, p)| = |Dp HM (x, p)| = I(x)f 2 |Dp FM (W (x, p))| 0 = f 2 kIk∞ |FM (W (x, p))Dp W (x, p)| . Hence, by (39) and (47), we get (H3) with mR (t) = Ct for any R ∈ (0, +∞). Validation of hypothesis (H4) For H(x, u, p) equal to H L and F equal to FL , we have H(x, u, p + λn(x)) ≤ 0 ⇒ F (W (x, p + λn(x))) ≤ e2R δf 2 for any u ∈ (−R, R), x ∈ Γ and p ∈ R2 . Recalling (37) and since F is invertible from [0, +∞) to [0, +∞), we get s   1 −1 e2R |p + λn(x)| ≤ F = CR . f2 δf 2 Therefore, we get (H4) for the Lambertian Hamiltonian H L . The previous argument works exactly in the same way for H(x, u, p) equal to either H P H or H BP since the corresponding functions FP H and FBP satisfy the same properties of FL (see Lemma 2.2(i)). This is not true for the Hamiltonian H ON defined in (15) since FON is bounded (see Lemma 2.2(ii)). Validation of hypothesis (H5) For H(x, u, p) equal to H L and F equal to FL , we have H(x, u, p − λn(x)) ≥ e−2R + δf 2 F (W (x, p − λn(x))) p = e−2R + δf 2 W (x, p − λn(x)) + 1 p ≥ e−2R + δf 2 f 2 |p − λn(x)|2 = e−2R + δf 2 f |p − λn(x)| 12 (50) that goes to +∞ for λ → +∞. Hence, we get (H5) for the Lambertian Hamiltonian (10). The previous argument applies also to the Phong Hamiltonian (22) since for W (x, p) ≥ 1 the behavior of H P H is as the one of H L . For the BP–model, considering the associated Hamiltonian H BP defined in (26) we have H BP (x, u, p − λn(x)) ≥ = ≥ = ≥ = e−2R + δf 2 F BP (W (x, p − λn(x))) c (W (x, p − λn(x)) + 1) 2 e−2R + δf 2 c−1 kD (W (x, p − λn(x)) + 1) 2 + kS c (W (x, p − λn(x)) + 1) 2 −2R 2 e + δf c−1 (kD + kS )(W (x, p − λn(x)) + 1) 2 p W (x, p − λn(x)) + 1 −2R 2 e + δf (kD + kS ) p f 2 |p − λn(x)|2 e−2R + δf 2 (kD + kS ) f |p − λn(x)| e−2R + δf 2 (kD + kS ) (51) that goes to +∞ for λ → +∞. Hence, we get (H5) for H BP . As a consequence of the previous Proposition 4.3 and Theorem 3.4, we get the following Corollary 4.4. We have (i) For H given by H L , H P H or H BP , the boundary value problem (45) with either Dirichlet, Neumann or state constraints boundary conditions in viscosity sense satisfies the Maximum Principle for discontinuous viscosity solutions and therefore it admits at most one viscosity solution. (ii) For H given by H ON , the boundary value problem (45) with Neumann boundary condition in viscosity sense satisfies the Maximum Principle for discontinuous viscosity solutions and therefore it admits at most one viscosity solution. We do not discuss here the existence of a viscosity solution to the boundary value problems (44) and (45). Since the SfS problem is an inverse problem in Computer Vision, we a priori know that, if we prescribe the correct boundary conditions, the solution u is given by the height of the surface to be reconstructed. Moreover, thanks to the Maximum Principle for discontinuous viscosity solutions, see Theorem 3.4, it is possible to get existence of a solution by proving the convergence of the various schemes considered in literature (see e.g. [7, 12]). Finally, the Maximum Principle is also a key ingredient to get existence of a solution via Perron’s method, once a subsolution and a supersolution of the problem have been constructed (see [9]). We conclude with some remarks about the different boundary conditions in Shape-fromShading theory. The choice of the boundary condition obviously depends on the chosen model and on available data. But, as observed in [22], the advantage of the state constraints boundary condition with respect to the Dirichlet and Neumann boundary conditions is that it does not require any a priori information about the solution at the boundary. Nevertheless, the approximation of this boundary condition can be difficult. In the next lemma we show that the solution of a Dirichlet problem with a large constant as boundary data and the solution of state constraints problem are upper bounded by the same constant M (independent of the solution). Therefore, in this case, by the very definition of boundary conditions in weak sense, the Dirichlet supersolution condition reduces to the state constraint one (46). It follows that to approximate a problem with state constraints boundary condition is sufficient to consider the corresponding problem with Dirichlet boundary datum given by a constant M sufficiently large. This request is not a strong constraint difficult to occur since, for example, if the input image contains an object of interest in front of a background, the condition is satisfied in a neighbourhood of the object where u(x) increases rapidly. 13 Lemma 4.1. For H(x, u, p) equal to either H L , H P H , or H BP , if the intensity image I in (48) verifies I(x) ≥ δ > 0, ∀x ∈ Ω, then the solution of the Dirichlet problem (44) with an appropriately large constant as boundary data and the solution of the problem (45) with state constrains conditions are upper bounded by the same constant (dependent only on the data of the problem). Proof. For H(x, u, p) equal to H L , let us define M L = − 21 ln(δf 2 ). Then, we have L H L (x, M L , 0) = −e−2M + I(x)f 2 ≥ 0 ∀x ∈ Ω. Therefore, the constant function u(x) ≡ M L is a supersolution of (44) in Ω. Moreover, ∀x ∈ ∂Ω and ∀ξ ∈ R2 , we have p L L H L (x, M L , ξ) = −e−2M + I(x)f 2 W (x, ξ) + 1 ≥ −e−2M + I(x)f 2 ≥ 0, (52) which implies that u(x) ≡ M L is a supersolution of (27) on the boundary and therefore, see (46), a supersolution of the state constraint boundary problem in Ø. Hence, by Theorem 3.4 and Corollary 4.4, the solutions of the Dirichlet and state constraints problems with Hamiltonian H L are upper bounded by the same constant M L . 2 We now prove the same result for the Phong model. Let us define M P H = − 12 ln( kDδf+kS ). We have H P H (x, M P H , 0) = = (I(x) − kA IA (x))f 2 kD + kS δf 2 (I(x) − kA IA (x))f 2 − + kD + kS kD + kS −e−2M PH + (53) (54) that is non-negative since I(x) = (I(x) − kA IA (x)) ≥ δ > 0, ∀x ∈ Ω. Therefore, the constant function M P H is a supersolution of (45) in Ω for the PH–model. Moreover, ∀x ∈ ∂Ω and ∀ξ ∈ R2 , we have for W (x, p) > 1, p PH H P H (x, M P H , ξ) = −e−2M + (I(x) − kA IA (x))f 2 W (x, ξ) + 1 ≥ −e−2M PH + (I(x) − kA IA (x))f 2 ≥ 0; for 0 ≤ W (x, p) ≤ 1, H P H (x, M P H , ξ) = −e−2M PH (W (x, ξ) + 1)α+1/2 kD (W (x, ξ) + 1)α + kS (W (x, ξ) + 1)1/2 (1 − W (x, ξ))α (I(x) − kA IA (x))f 2 + ≥ 0. kD + kS + (I(x) − kA IA (x))f 2 ≥ −e−2M PH So, the constant function M P H for the PH–model is a supersolution of (45) in Ω for both the Dirichlet problem with boundary datum M P H and the state constraints one. For the BP–model, we can consider the same constant used for the PH–model and repeat the same steps to prove the Lemma. 5 Conclusions In this paper we have shown the well-posedness of several non-Lambertian models in an analytical way thanks to the introduction of an attenuation term, under a perspective projection with a unique light source located at the optical center of the camera. We have derived the Hamilton-Jacobi equations associated to each reflectance model and then we have obtained a 14 unified formulation showing useful properties for the general Hamiltonian. We have clarified the differences between the various models from the point of view of the problem with Dirichlet boundary conditions (BC), Neumann BC and state constraints BC. This analysis explains (and corrects) in details the empirical deductions coming from previous works on non-Lambertian models (see e.g. the conclusion reported in [1] regarding the concave/convex ambiguity for the ON-model), focusing the attention on the important role of the boundary conditions and providing also an useful theoretical tool for the numerical point of view thanks to the last Lemma 4.1. It is also worth noting that, from a theoretical point of view, contrary to previous works (see e.g. [11, 18, 15, 10]), we do not require any additional information for the models and we do not need to add regularization terms to the PDEs for ensuring the well-posedness of the SfS problem. From a numerical point of view, the unified approach allows to consider a unique and general numerical scheme that incorporates all the different SfS models described in the paper (and possibly other non-Lambertian models not considered here), thereby providing an easy way to compare the performances of the different models by only changing some parameters. A Proofs of lemmas in Subsection 2.5 Proof of Lemma 2.1. We have f 2 |p|2 ≤ W (x, p) = |p|2 W (x, p ) ≤ |p|2 max W (x, q) |p| |q|=1 which proves (37). Moreover, i |x|2 + f 2 h 2 2 2 2x + (p · x) + f |p| f2 f2 2 2 |x| + f Dp W (x, p) = 2(f 2 + |x|2 )p + 2 (x · p)x f2 Dx W (x, p) = 2(p · x)p from which we get (38) and (39).  Proof of Lemma 2.2. (i) For the function FL we have lim FL (r) = lim r→+∞ r→+∞ √ 1 r + 1 = +∞, FL0 = √ 2 r+1 for r ∈ [0, +∞), and the corresponding properties in the statement easily follow. Regarding the PH–model, we analyze the function FP H (r) only for r ∈ [0, 1) (otherwise it is FL (r) equal to ). For r ∈ [0, 1), we have kD 1 · [kD (r + 1)α + kS (r + 1)1/2 (1 − r)α ]2 h 1   ( + α) kD (r + 1)2α−1/2 + kS (r + 1)α (1 − r)α − αkD (r + 1)2α−1/2 2 1 i −kS (r + 1)α (1 − r)α − α(r + 1)α+1 (1 − r)α−1 2 i h1 1 2α−1/2 α α−1 = k (r + 1) + 2αk (r + 1) (1 − r) . · D S [kD (r + 1)α + kS (r + 1)1/2 (1 − r)α ]2 2 FP0 H (r) = Hence, h1 i 1 2α−1/2 α k (r + 1) + 2αk (r + 1) D S 2 (r + 1)2α 2 kD 1 2αkS √ = + 2 k (r + 1)α 2kD r + 1 D  1  2αkS 1 √ ≤ + 2 . 2kD kD r+1 0 < FP0 H (r) ≤ 15 Moreover, lim FP0 H (r) = lim+ FP0 H (r) = r→1− r→1 1 0 F (1), kD L hence FP H ∈ C 1 . By the corresponding properties for FL , it follows that FP H is strictly increasing and the following bounds hold 0 < FP0 H (r) ≤ C √ 1 r+1 for r ∈ [0, +∞) and ∀α ≥ 1. 1 For the BP–model, if c = 1 in (36), then FBP (r) = kD +k FL (r), hence we can refer to FL . In S the case c > 1, the function FBP is still smooth and positive. Moreover, lim FBP (r) = +∞ r→∞ and 0 FBP (r) = ≤ c 2 (r c + 1) 2 −1 (kD (r + 1) c−1 2 c + kS ) − (r + 1) 2 ( c−1 2 kD (r + 1) (kD (r + 1) c 2 (r + 1) c 2 −1 (kD (r + 1) (kD (r + 1) c−1 2 c−1 2 + kS ) + kS )2 c−1 2 = + kS c−1 2 −1 ) )2 c 2 (r c + 1) 2 −1 kD (r + 1) c−1 2 + kS . Hence, the following bounds hold 0 0 < FBP (r) ≤ C √ 1 r+1 for r ∈ [0, +∞) and ∀c ≥ 1. This shows that also the function FBP is strictly increasing, that concludes the proof of (i). (ii) For the function FON we have √ A r + 1 − 2B 0 √ for r ∈ [0, +∞). FON (r) = 2[A r + 1 + Br]2 Hence, FON is strictly increasing if A/2 > B and the property (43) in the statement easily follows. Moreover, since limr→+∞ FON (r) = 1/B, it is bounded.  References [1] A.H. Ahmed and A.A. Farag. A New Formulation for Shape from Shading for NonLambertian Surfaces. In Proc. IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 1817–1824, 2006. [2] A.H. Ahmed and A.A. Farag. Shape from Shading Under Various Imaging Conditions. In IEEE International Conference on Computer Vision and Pattern Recognition CVPR’07, pages X1–X8, Minneapolis, MN, June 2007. [3] M. Bardi and I. Capuzzo-Dolcetta. Optimal Control and Viscosity Solutions of HamiltonJacobi-Bellman Equations. Modern Birkhäuser Classics. Birkhäuser, Boston, 2008. [4] G. Barles. Solutions de viscosité des équations de Hamilton-Jacobi, volume 17 of Mathématiques & Applications (Berlin) [Mathematics & Applications]. Springer-Verlag, Paris, 1994. [5] G. Barles and P.E. Souganidis. Convergence of approximation schemes for fully nonlinear second order equations. Asymptotic Anal., 4(3):271–283, 1991. [6] J.F. Blinn. Models of Light Reflection for Computer Synthesized Pictures. Computer Graphics, 11(2):192–198, 1977. 16 [7] M. Breuß, E. Cristiani, J.-D. Durou, M. Falcone, and O. Vogel. Perspective Shape from Shading: Ambiguity Analysis and Numerical Approximations. SIAM J. Imaging Sciences, 5(1):311–342, 2012. [8] A.R. Bruss. The Image Irradiance Equation: Its Solution and Application. Technical report ai-tr-623, Massachusetts Institute of Technology, June 1981. [9] M.G. Crandall, H. Ishii, and P.L. Lions. User’s Guide to Viscosity Solutions of Second Order Partial Differential Equations. Bull. Amer. Math. Soc. (N.S.), 27(1):1–67, 1992. [10] A. Crouzil, X. Descombes, and J. D. Durou. A Multiresolution Approach for Shape from Shading Coupling Deterministic and Stochastic Optimization. IEEE Transactions on Pattern Analysis and Machine Intelligence, 25(11):1416–1421, Nov 2003. [11] Paul Dupuis and John Oliensis. An Optimal Control Formulation and Related Numerical Methods for a Problem in Shape Reconstruction. Ann. Appl. Probab., 4(2):287–346, 05 1994. [12] J.-D. Durou, M. Falcone, and M. Sagona. Numerical methods for shape-from-shading: A new survey with benchmarks. Computer Vision and Image Understanding, 109(1):22–43, 2008. [13] M. Falcone and R. Ferretti. Semi-Lagrangian Approximation Schemes for Linear and Hamilton-Jacobi Equations. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2014. [14] B.K.P. Horn and M.J. Brooks. The Variational Approach to Shape from Shading. Computer Vision, Graphics and Image Processing, 33(2):174–208, 1986. [15] B.K.P. Horn and M.J. Brooks. Shape from Shading. MIT Press, Cambridge, MA, USA, 1989. [16] Y.-C. Ju, S. Tozza, M. Breuß, A. Bruhn, and A. Kleefeld. Generalised Perspective Shape from Shading with Oren-Nayar Reflectance. In Proceedings of the 24th British Machine Vision Conference (BMVC 2013), pages 42.1–42.11, Bristol, UK, 2013. BMVA Press. [17] T. Okatani and K. Deguchi. Shape Reconstruction from an Endoscope Image by Shape from Shading Technique for a Point Light Source at the Projection Center. Computer Vision and Image Understanding, 66(2):119–131, 1997. [18] J. Oliensis and P. Dupuis. A Global Algorithm for Shape from Shading. In Computer Vision, 1993. Proceedings., Fourth International Conference on, pages 692–701. IEEE International Conference on Computer Vision, May 1993. [19] M. Oren and S.K. Nayar. Generalization of the Lambertian Model and Implications for Machine Vision. International Journal of Computer Vision, 14(3):227–251, 1995. [20] M. Oren and S.K. Nayar. Generalization of Lambert’s Reflectance Model. In SIGGRAPH ’94 Proceedings of the 21st annual conference on Computer graphics and interactive techniques, pages 239–246, Orlando, FL, USA, July 1994. [21] B.T. Phong. Illumination for Computer Generated Pictures. Communications of the ACM, 18(6):311–317, 1975. [22] E. Prados and O. Faugeras. Shape From Shading: a well-posed problem? In International Conference on Computer Vision and Pattern Recognition CVPR05, volume 2, pages 870 – 877, San Diego, CA, USA, June 2005. [23] E. Prados, O. Faugeras, and F. Camilli. Shape from Shading: a well-posed problem? Technical Report 5297, INRIA, August 2004. [24] E. Rouy and A. Tourin. A Viscosity Solutions Approach to Shape-From-Shading. SIAM J. Numer. Anal., 29(3):867–884, 1992. [25] S. Tozza and M. Falcone. Analysis and Approximation of Some Shape–from-Shading Models for Non-Lambertian Surfaces. Journal of Mathematical Imaging and Vision, 55(2):153–178, 2016. [26] O. Vogel, M. Breuß, and J. Weickert. Perspective Shape from Shading with Non-Lambertian Reflectance. In Gerhard Rigoll, editor, Pattern Recognition, volume 5096 of Lecture Notes in Computer Science, pages 517–526. Springer Berlin Heidelberg, 2008. 17 [27] R. Zhang, P.-S. Tsai, J.E. Cryer, and M. Shah. Shape from Shading: A Survey. IEEE Transactions on Pattern Analysis and Machine Intelligence, 21(8):690–706, 1999. 18