Perspective Shape-from-Shading Problem: A Unified Convergence Result for Several Non-Lambertian Models

Shape-from-Shading represents the problem of computing the three-dimensional shape of a surface given a single gray-value image of it as input. In a recent paper, we showed that the introduction of an attenuation factor in the brightness equations related to various perspective Shape-from-Shading models allows us to ensure the well-posedness of the corresponding differential problems. Here, we propose a unified convergence result of a numerical scheme for several non-Lambertian reflectance models. This result is interesting since it can be easily extended to other non-Lambertian models in a unified and, therefore, powerful framework.


Introduction
The Shape-from-Shading (SfS) problem consists of computing the three-dimensional shape of an object starting from a single gray-level image of it. This problem aroused the interest of some opticians in the 1950s-1960s [1,2], and then the problem was formulated by B.K.P. Horn and his collaborators of MIT as a first-order non-linear partial differential equation (PDE) in the 1970s-1980s [3][4][5]. The goal was to enable the 3D surface represented in the input image to solve a Hamilton-Jacobi (HJ) equation or a variational problem. This problem gave rise to an expansion in the field of mathematics, and some researchers tried to prove the well-posedness of it in the framework of weak solutions. The first works of Lions, Rouy and Tourin in the early 1990s [6,7] inserted the SfS problem into the context of the viscosity solutions frameworks, hence in a much more theoretical area.
In general, the SfS problem is described by the following image irradiance equation [8]: In this equation, the gray-level intensity I(x) of the given input image at the point x := (x, y) is put in relation to the reflectance function R(N(x)), which represents the value of the light reflection on the surface as a function of its orientation, i.e., of the unit normal N(x). Different expressions of the function R produce different reflectance models.
There are three main ingredients for formulating the SfS problem: how we model the light, the camera and the reflectance. The most common setup found in the literature considers the Lambertian reflectance model (for which the intensity I(x) depends only on the incident angle between the direction of the light source ω and the direction of the normal unit vector, without considering the viewer), under an orthographic projection (no perspective deformations are considered) with a single light source located far from the object to be reconstructed (located at infinity in the direction of the light versor ω). In this case, the irradiance Equation (1) can be written as the following HJ equation: where u is the unknown surface height and γ D is the diffuse albedo, which indicates the reflective power of a diffuse surface. In general, non-linear PDEs such as (2) do not admit smooth solutions. Hence, we need to look for a proper notion of a weak solution, for example, in the context of the viscosity ones (see [9][10][11]). Unfortunately, under orthographic projection, there is no uniqueness of solution, not only in the classical sense but also in the weak sense, due to the well-known concave/convex ambiguity related to singular points, i.e., points where the intensity of the image is at the maximum and the light direction is parallel to the surface normal, generating countless viscosity solutions to Equation (2). In order to overcome this ambiguity, several attempts have been made by the scientific community in different research directions. Considering a priori information on the object shape model and/or an integrability condition of the height field in order to constrain the ambiguities could be a possibility. From the mathematical point of view, adding information considering more than one input image, taken from the same point of view but with different light sources (photometric-stereo technique [12,13]) or under the same light source but with different viewpoints (stereovision; see [14]), allows us to achieve well-posedness. In order to stay within the scope of a single input image, a perspective Lambertian model was proposed in [15][16][17], with a light source located far from the surface to reconstruct, i.e., located at infinity. In [16], and then in the extended journal version [17] by the same authors, it was observed that the perspective projection is a more realistic projection than the orthographic one, and it improves the performance with respect to considering the orthographic projection, regardless of the specific numerical algorithm employed (in [17], a comparison with [15] and three orthographic SfS methods is also reported, in order to show the better performance of the perspective models). In [18,19], a different setup for the perspective model was proposed, with a pinhole camera and light source located at the optical center. All these works are limited to the assumption of a Lambertian reflectance model, which is known to be not always suitable for describing real-world surfaces. Models for non-Lambertian surfaces (e.g., [20][21][22]) have been proposed in chronologically later works [23][24][25][26]. All these models, when formulated in terms of a differential equation, do not resolve the concave/convex ambiguity (see [27] for an analysis of them). In order to achieve the uniqueness of solution starting from a single input image, in [18,28], the authors introduced an attenuation factor in the Lambertian brightness equation under perspective projection. This factor takes into account the distance between the surface and the light source and, thanks to this attenuation term, the associated HJ equation admits a unique viscosity solution. Starting from the idea in [18] related to the Lambertian model, in [29], the authors showed the well-posedness of several non-Lambertian models in a unified formulation, by validating the assumptions of the Maximum Principle for discontinuous viscosity solutions (see [9]). Many works have focused on the approximation of the different brightness equations derived from the various models considered for solving the SfS problem (the reader can refer to the surveys [30,31] for more details), but only very few of them treat the convergence of the corresponding methods, and the results are often of experimental type. In this paper, we show the convergence of a numerical scheme for the perspective SfS problem associated with different non-Lambertian models, based on the method proposed in [32] and the theoretical results contained in [29]. To the best of our knowledge, this is the first paper in which a theoretical convergent result of a numerical scheme related to different non-Lambertian models in a unified framework is provided.
The paper is organized as follows: in Section 2, we introduce the setup and the notation adopted, briefly recalling the reflectance models we considered and how it is possible to write them in a unified formulation. Then, we recall the definition and basic properties of viscosity solutions and the theoretical results shown in [29]. Section 3 represents the core of this work, in which we state and prove the convergence result. In Section 4, numerical experiments confirm the effectiveness of the proposed scheme. The paper ends with final comments and conclusions reported in Section 5.

Preliminary Materials and Models
In this section, we introduce the notation and the setup we adopted, before recalling the four reflectance models considered. Finally, we introduce the definition of viscosity solutions and their basic properties.

Setup and Notation
For all four models, for the setup, we used the perspective camera projection with one light source located at the optical center of the camera and an attenuation term of the illumination due to the distance between the surface and the light source.
Let Ω ⊂ R 2 be open and bounded, which represents the image domain. As performed in [18,23,25,28], the scene can be represented by a surface S : Ω → R 3 defined as the following function: where u(x) is the unknown height of the surface and f > 0 is the focal length of the camera. The unit vector in the direction of the light source is defined as follows: The unit normal vector at the point x of S is defined as and

Reflectance Models
We now briefly recall the reflectance models and their corresponding HJ equations, describing them in a unified formulation (see [29] for more details).

Lambertian Model
The Lambertian model (L-model) is the most common and simple reflectance model ( [4,5,18,19]), which does not consider the viewer position. The brightness equation for this model, considering an attenuation term, is as follows: where I(x) denotes the intensity of the image, γ D (x) is the diffuse albedo, θ i is the incident angle, i.e., the angle between the light source direction ω(S (x)) and the unit normal to the surface N(x), and 1/r 2 is the attenuation term of the illumination due to the distance between the surface and the light source [23]. Under the setup introduced in Section 2.1, supposing γ D (x) = 1, i.e., all incident light is reflected, and considering the attenuation term 1/r 2 , with r = f u(x), the following HJ equation can be associated with (7): where By defining we can rewrite Equation (8) as follows: Supposing that the surface S is in front of the optical center, i.e., it is visible, then u is a strictly positive function. Hence, we can consider the change of variable v(x) = ln u(x), arriving at the following HJ equation in the new variable v (cf. [18,28]) to which we associate the Hamiltonian

Oren-Nayar Model
Considering real-world surfaces, e.g., plaster, sand, or concrete, the Lambertian model is an inadequate approximation of the diffuse component, mainly because it does not take into account the roughness of the surface. In the 1990s, Michael Oren and Shree K. Nayar proposed a different diffuse model, suitable for rough surfaces, which are modelled as a set of facets with different slopes, each of them Lambertian in reflectance [22,33]. This model takes into account complex physical phenomena, such as masking, shadowing and interreflections between points on the surface facets, and can be viewed as a generalization of the Lambertian one. The brightness equation for the Oren-Nayar model (ON-model) is as follows: where θ i indicates the incident angle, θ r denotes the angle between the observer V and the normal N directions, α and β are defined as α = max{θ i , θ r } and β = min{θ i , θ r }, respectively. ϕ i and ϕ r are, respectively, the angles between the projection of the light source direction ω, or the projection of the viewer direction V, and the x axis onto the (x, y)plane. The constants A and B are non-negative quantities and depend on the roughness parameter σ as follows: As a consequence of the chosen setup illustrated in Section 2.1, θ i = θ r = α = β, and we can simply call it θ. Hence, adding the attenuation factor 1/r 2 as proposed in [23], the Oren-Nayar brightness equation is as follows: Note that when σ = 0, then A = 1 and B = 0, so the ON-model goes back to the L-model. Since r = f u(x), and cos(θ) = N(x) · ω(S (x)), from which sin 2 (θ) = 1 − (N(x) · ω(S (x))) 2 , then the brightness Equation (14) becomes: where Q(x) is defined as in (9). By using the function W defined in (10), and the change of variable v(x) = ln u(x) (which is possible since u(x) > 0 ∀x ∈ Ω), we obtain the following HJ equation in v: and we associate with it the following corresponding Hamiltonian . (17) Note that the function W(x, p), defined in (10), is the same function used in the L-model, which is a useful tool for the unified approach of different reflectance models.

Phong Model
The Phong model (PH-model) [20] belongs to the category of specular reflectance models. The law of specular reflection states that the angle of reflection of a ray equals the angle of incidence, and that the incident direction, the surface normal and the reflected direction are coplanar. The PH-model considers a specular term in the definition of the function I(x): where I A (x) is the ambient light component and k A the percentages of this component, the diffuse component is defined as in the L-model, with k D indicating the percentages of this component, and the third addend represents the specular light component, with k S indicating the percentages of it. This last term is described as a power of the cosine of the angle θ s between the unit vectors V and R(x), with R indicating the reflection of the light ω on the surface; γ S (x) denotes the specular albedo, and α represents the characteristics of specular reflection of a material. We assume k A + k d + k S = 1, and that α is an integer greater or equal to 1 (in [20], Phong assumes α ∈ [1,10]). Under the chosen setup, the direction of the light source and that of the viewer are the same; therefore, θ s = 2θ i (see [25]). By setting the diffuse and specular albedo to 1 and adding the light attenuation term, we can arrive at the following non-linear PDE associated with the brightness Equation (18): where Q(x) is defined as in (9), |n(x)| as in (6), (10), we can obtain the following HJ equation in v: −e −2v(x) A different equation was derived in [25], defining an alternative function W, which, however, is not useful for our purpose of a unified formulation.
As stated in [25], since the term 1+W(x,p) represents the cosine of the specular term, in the PH-model, this cosine is replaced by zero if cos θ s = 1−W(x,p) 1+W(x,p) < 0. Hence, we can associate to (20) the following Hamiltonian: Note that the Phong model reduces to the Lambertian one (up to a constant) if k S = 0, i.e., the specular component vanishes.

Blinn-Phong Model
In 1977, Blinn [21] proposed a modification of the Phong model by introducing an intermediate vector H, which bisects the angle between the unit vectors V and ω. If the surface is a perfect mirror, the light reaches the viewer only if the surface normal N is pointed halfway between the viewer direction V and the light source direction ω. The direction of maximum highlight is denoted by H = ω+V |ω+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 H and N 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 power (in [21], it is stated that this power is typically 50 or 60).
For the Blinn-Phong model (BP-model), the brightness equation is as follows: where δ is the angle between H and the unit normal N, and c measures the shininess of the surface (we assume c ≥ 1). Using the setup introduced in Section 2.1, the directions of the light source and the viewer are the same. As already carried out for the PH-model, we set γ D (x) = γ S (x) = 1 and we add the light attenuation term as performed for the previous models. In this way, Equation (22) becomes the following: ) and cos δ = N(x) · ω(S (x)) due to the chosen setup. Using the definition of the unit normal N(x), thanks to (5) and (6), the function W(x, ∇v(x)) as defined in (10), Q(x) as defined in (9), and the same change in variable adopted before, i.e., v(x) = ln u(x), we obtain the following HJ equation in the new variable v (see [29] for more details): with which we associate the Hamiltonian Remark 1. As already observed for the PH-model, if k S = 0, the BP-model also reduces to the Lambertian one, up to a constant.

Reflectance Models in a Unified Formulation
The four different reflectance models presented in Section 2.2 lead to a HJ equation of the following form: where v(x) is related to the original unknown u(x), which represents the surface height, thanks to the change of variable v(x) = ln u(x). The Hamiltonian H has a peculiar explicit formulation for each of the four models, given by (12), (17), (21) and (24), respectively. Nevertheless, looking at the four Hamiltonians, we can note that (25) can be rewritten separating a term depending only on v(x) and a Hamiltonian which depends only on (x, ∇v(x)) as follows: where M indicates the acronym of the four models (M = L, ON, PH, BP). Since the function W(x, p) adopted is the same for all four cases, then H M (x, p) corresponds for each one to the following: Lambertian model where Oren-Nayar model where Phong model where where In [29], some useful bounds on the function W were proved. In particular, the authors showed the following: where the constant C is independent of (x, p) ∈ Ω × R 2 . These bounds help to prove our convergence result. Some properties of the four functions F M can also be proved (see Lemma 4 in [29]).

Viscosity Solutions of Hamilton-Jacobi Equations
Let us recall in this section the definition and basic properties of viscosity solutions in order to study the general Hamiltonian (25) in the case of our reflectance models (the reader can refer to [9,11] for more details on this theory). As already stated in the introduction, a HJ equation such as (25), in general, does not admit a classical solution, and we need to look for a proper notion of a weak solution, for example, in the context of viscosity solutions. However, if a classical solution to the HJ equation exists, it coincides with the viscosity one. Tipically, viscosity solutions are Lipschitz continuous functions, which can be obtained via the method of vanishing viscosity, i.e., as a limit of regular solutions to second-order problems.
We can now give the definition of a viscosity subsolution, supersolution and solution.
Finally, v is a viscosity solution of (25) if it is both a viscosity super-and subsolution.
To prove the convergence of a numerical scheme using the Barles-Souganidis approach [32], we need the Maximum Principle for discontinuous viscosity solution, which requires the formulation in a weak (viscosity) sense of the boundary conditions for the PDE to be approximated. Of course, this requirement is in addition to the consistency, monotonicity and stability properties of the numerical scheme. Hence, motivated by this remark, we need to introduce the notion of boundary conditions in a weak (viscosity) sense, which means we have to incorporate the boundary condition into the definition. To gain a better understanding, instead of considering, for example, the Dirichlet boundary condition in a pointwise sense associated with the HJ Equation (25), i.e., v(x) = g(x) ∀x ∈ ∂Ω, we introduce the operator B : ∂Ω × R × R 2 → R, which represents the boundary condition, and in the case of the Dirichlet boundary condition, it is chosen as B(x, s, p) = s − g(x). Let us consider a more general boundary value problem of the following form: with the operator B : ∂Ω × R × R 2 → R. We now introduce the notion of (discontinuous) viscosity solution for (35).

Definition 2.
A viscosity supersolution of (35) is a LSC function v in Ω s.t. for any ψ ∈ C 1 (Ω), if x 0 is a local minimum point of v − ψ, then A viscosity subsolution of (35) is a USC function v in Ω s.t. for any ψ ∈ C 1 (Ω), if x 0 is a local maximum point of v − ψ, then Finally, v is a viscosity solution of (35) if it is both a viscosity super-and subsolution.
We are now ready to state the Maximum Principle for discontinuous viscosity solutions. We first consider the case of a Dirichlet boundary condition (see [9]).

Theorem 1.
Let Ω ⊂ R 2 be bounded with ∂Ω of class W 2,∞ , g : ∂Ω → R a continuous function and H : Ω × R × R 2 → R satisfying the following two conditions: Moreover, let us assume that there exists a neighborhood Γ of ∂Ω in R 2 , such that Let v be a USC function in Ω, subsolution of (35) (u is a LSC function in Ω, supersolution of (35), respectively) with B(x, s, p) = s − g(x). Then, v ≤ u in Ω.
This theorem implies the uniqueness of the viscosity solution in Ω, i.e., if v 1 , v 2 ∈ C(Ω) are two viscosity solutions of (35), then v 1 = v 2 in Ω. If we want to consider different boundary conditions and validate the Maximum Principle, it is sufficient to choose the following: and assuming that J (x) is a Lipschitz continuous function in Ω and there exists a neighborhood Γ of ∂Ω, such that J (x) ≥ δ > 0 in Γ, in [29], the authors validated the assumptions for the Maximum Principle related to the four Hamiltonians, H L , H ON , H PH , and H BP defined in (12), (17), (21) and (24), respectively, arriving to prove the well-posedness of the perspective SfS problems.

Remark 2.
Without the introduction of the attenuation term and/or under a different setup, e.g., in the context of orthographic projection, the assumption (36) would not be satisfied; hence, uniqueness would not be guaranteed.

Convergence Result
Following [32], we now establish necessary notions for the convergence of numerical schemes in the viscosity sense: consistency, monotonicity and stability. Afterwards, we bring in the convergence theorem [32], with which we can prove the convergence of our scheme.
To approximate (25), we consider a scheme of the following form: where S : , ρ is a discretization parameter and B(Ω) denotes the space of bounded functions defined on Ω. (42) is assumed to be consistent, provided that both the following conditions are satisfied:

Definition 3 (Consistency). The scheme defined by
For all x 0 ∈ Ω and ψ ∈ C 1 (Ω) and lim sup where H * and H * denote viscosity sub-and supersolution, respectively, defined in the previous section. (42) is assumed to be monotone, provided that the following condition is fulfilled S(ρ, x, t, u) ≤ S(ρ, x, t, v)

Definition 4 (Monotonicity). The scheme S in
if u ≥ v for all ρ ≥ 0, x ∈ Ω, t ∈ R and u, v ∈ B(Ω). (42) is called stable if for all ρ > 0, there exists a solution u ρ ∈ B Ω of (42) with a bound that is independent of ρ.

Definition 5 (Stability). The scheme S in
Based on these definitions, the convergence theorem in the viscosity framework is elaborated by the following fundamental result (cf. [32]): (42) is consistent, monotone and stable, then the solution u ρ of S converges locally uniformly to the unique viscosity solution u of the original problem (25) as ρ → 0.

Numerical Discretization
For the sake of simplicity, we stick to the 1-D version of (26) since this suffices to perceive the methodology. The 1-D model of (26) reads as follows: where M is the acronym of the models (M = L, ON, PH, BP). For each model, the corresponding H M (x, v x (x)) uses the 1D version of (9), (10) and the various functions F M .
Let ∆x > 0 be the spatial step size, and we denote by N := N(∆x) the number of grid points x i , i = 1, . . . , N. We further denote the approximate value of v at the i-th grid point x i by w i and define ψ i (w) as the following: where w = (w 1 , . . . , w N ), cf. [34]. In addition, for consistency, the approximate spatial derivative is given by the following: In order to discretize (46), we introduce the method of artificial time, and we make use of the Euler forward in time scheme in addition to (48), obtaining at each grid point the following approximation scheme: which we can rewrite to obtain the following update rule: For simplicity, here, we also assume that the iteration level used for the discrete representations of v x is always the actual level n. This also holds for the source term, i.e., e −2v ≈ e −2w n i .

Proof of the Convergence Result
We now state the main result of this paper: Theorem 3. Let us assume that the scheme given in (50) satisfies the assumptions of the Theorem 2 and the HJ equation in (46) admits a strong uniqueness property. Then, the solution of (50) converges locally uniformly to the unique continuous viscosity solution of (46).

Proof.
We have to validate the assumptions of strong uniqueness, consistency, monotonicity and stability. Regarding the strong uniqueness for (46), Propositions 9 and 11 in [29] already validated the required conditions. In addition, the scheme (50) is consistent with the HJ Equation (46) since this is directly followed by the construction. Moreover, by [35], the stability holds when the scheme is monotone. Hence, the verification of the monotonicity remains. To this end, by [35], it is sufficient to show that G in (50) is a non-decreasing function, which holds if and only if: Due to analogy between case (i) and (iii), we proceed in the following order: case (i) → (iii) → (ii).

Case (i).
In view of (47) and (48), there is a contribution to the following case (i) only when the third argument in (47) is taken. Taking a partial derivative of (50) with respect to w i−1 yields the following: Note that −∆t(I f 2 ) or −∆t(I − I a k a ) f 2 are both negative quantities. Therefore, in At this point, we need to distinguish the four cases. Note that for construction and for Lemma 3 of [29], the function W(x, v x ) ≥ 0. In the following, we omit the variable dependencies for brevity.

Oren-Nayar case (F ON
Since the denominator is always positive, we focus only on the numerator: Since the denominator is always greater than zero, we continue focusing only on the numerator.
In the case of W ≥ 1, we fall back on the Lambertian case, so it is already proven. In the case 0 ≤ W < 1 Since the denominator is always greater than zero, we continue focusing only on the numerator.
≥ 0 also for the Phong case.

Case (iii).
For the case (iii), there is a contribution only when the second argument in (47) is chosen. Taking partial derivative of (50) with respect to w i+1 yields the following: As for case (i), note that −∆t(I f 2 ) or −∆t(I − I a k a ) f 2 are both negative quantities.
Therefore, in order to obtain ∂G ∂w i+1 ≥ 0, it is necessary to show that ∂F M ∂w i+1 ≤ 0.
We note that ∂W ∂w i+1 Hence, the quantity ∂W ∂w i+1 is less than zero in this case, as it was ∂W ∂w i−1 in the case (i).
As a consequence of this remark, the proof for this case (iii) is analogous to the previous case (i) for all four functions F, resulting in ∂G ∂w i+1 ≥ 0 for each model under the same assumptions.

Case (ii).
Taking partial derivatives of (50) with respect to w i results in the following: At this point, we need to distinguish the four models.
recalling thatŵ x,i := max(0, −D + w, D − w). By definition ofŵ x,i , we note that it is always, in any case, a non-negative quantity and Since a violation of case (ii) primarily occurs due to the quantity f 2 + x 2 ŵ x,i in the last term of (63), we incorporate a worst case estimate as follows: where the maximum is taken over all possible cases in (48). It follows Hence, we can formulate a sufficient condition for monotonicity on the time step size as follows: As a consequence, for the Lambertian case, the scheme satisfies monotonicity and, therefore, stability as long as (67) is valid.
Recalling the definitionŵ x,i := max(0, −D + w, D − w) and the estimate (65) used in the previous Lambertian model case, we arrive at the following condition: Hence, for the Oren-Nayar case, the scheme satisfies monotonicity and, therefore, stability as long as (69) is valid. We start focusing on ∂F BP ∂w i : Since the last term is positive, recalling the estimate (65) whereŵ x,i := max (0, −D + w, D − w), we can finally write the following: Hence, which guarantees monotonicity and, therefore, stability for the Blinn-Phong case.

Phong case
In the case of W ≥ 1, we fall back on the Lambertian case, which we have already seen.
Hence, we obtain Since the last term is (. . . ) ≥ 0, recalling the definition ofŵ x,i and the estimate (65), we can obtain the following: Hence, the monotonicity and, therefore, stability is also guaranteed for the Phong case.

Numerical Simulations
In this section, we report some one-dimensional synthetic numerical experiments, which confirm the convergence of the proposed scheme to a viscosity solution. All the numerical tests were implemented in Matlab, by using a Notebook HP EliteBook 830 G6 Intel Core i5-8265U with a speed of 1.60 GHz and 16 GB of RAM.

Test 1.
The unknown function u we want to reconstruct is defined by the following equation: Since it is a synthetic case, the input image is different for each model, generated by the proper brightness equation. This is clearly visible in Figure 1, where the true and the approximated images for each model are reported. The approximated four image functions I are computed a posteriori using the u computed by the scheme and central finite difference to obtain the gradient of u.
In Figure 2, we can see the exact and the approximated scenes defined according to Equation (3). For all the models, the approximation is really good.
Finally, in Figure 3, we report the plots of the exact and approximated heights of the surfaces related to each of the four reflectance models. Qualitatively, in all four cases, the solution computed by the scheme seems to be a good approximation of the surface height, confirming that the numerical solution converges to the correct (and unique) viscosity solution. In fact, looking at the pictures reported in Figure 3, it is difficult to distinguish the exact solution from the approximate one, except for in the PH-model, where it differs going towards the boundary. In order to also perform a quantitative analysis of these results, we compute the Root Mean Square Errors (RMSEs) related to the solutions of the four cases for each point, as shown in Figure 4 and defined as follows: Analyzing those figures, we can state that the RMSEs are of order 10 −5 for the L-model and BP-model, whereas they are are of order 10 −4 for the ON-model and PH-model.
Looking at Figures 5 and 6, we can see that an ambiguity, which remembers a concave/convex ambiguity typical of an orthographic setup, is present on the right side of the scene, even if the image is perfectly reconstituted, as shown, for example, in Figure 7 for the Phong and Blinn-Phong models. Looking at Figure 8, we can note that for both models the error is close to zero, except on the right side of the pictures, in correspondence to the ambiguity, where there is an error of about 0.01.
The ambiguity of this particular case is due to the model, not to the specific numerical scheme chosen to solve the problem. The existence of different surfaces which share the same image does not contradict the uniqueness result related to the viscosity solution, since they can correspond to different boundary conditions, or to the same boundary conditions imposed in a different domain of definition. This is just to stress that the uniqueness of the viscosity solution does not completely solve the problem of model ambiguity, since we could be interested in reconstructing a surface described by a weak solution, which is different from the viscosity solution. We refer the reader to [36] for more details on the ambiguities of the model. The mentioned paper refers only to the Lambertian case, but the analysis performed, even if only on a single model, is enough to understand what happens.

Conclusions
In this paper, we have shown a convergence result for the recent unified formulation of the perspective SfS models proposed in [29], in which the authors considered an attenuation term in order to achieve the well-posedness of the problem in the context of viscosity solutions. Despite the peculiarities of the different reflectance models considered, by formulating them in a unified framework described by a general Hamiltonian, which has attractive properties, we proved that a monotone, consistent and stable explicit scheme of upwind type converges to a viscosity solution, under uniqueness properties of the Hamiltonian. This result is interesting and powerful since can be easily extended to other non-Lambertian models in a unified framework. The numerical experiments shown confirm the convergence of the proposed scheme. We believe that this work could help to develop convergent numerical methods of refined SfS processes for real-world scenarios in future research.
Funding: This research received no external funding.

Informed Consent Statement: Not applicable.
Acknowledgments: The author thanks the anonymous referees for their useful suggestions, which allowed to improve the original version of the paper. The author is a member of the INdAM Research Group GNCS.

Conflicts of Interest:
The author declares no conflicts of interest.