The 3D reconstruction of an object starting from one or more images is a very interesting inverse problem with many applications. In fact, this problem appears in various fields which range from the digitization of curved documents  to the reconstruction of archaeological artifacts . More recently, other applications have been considered in astronomy to obtain a characterization of properties of planets and other astronomical entities [66, 27, 41]
and in security where the same problem has been applied to the facial recognition of individuals.
In real applications, several light sources can appear in the environment and the object surfaces represented in the scene can have different reflection properties because they are made by different materials, so it would be hard to imagine a scene which can satisfy the classical assumptions of the 3D reconstruction models. In particular, the typical Lambertian assumption often used in the literature has to be weakened. Moreover, despite the fact that the formulation of the Shape-from-Shading problem is rather simple for a single light source and under Lambertian assumptions, its solution is hard and requires rather technical mathematical tools as the use of weak solutions to nonlinear partial differential equations (PDEs). From the numerical point of view the accurate approximation of non regular solutions to these nonlinear PDEs is still a challenging problem. In this paper we want to make a step forward in the direction of a mathematical formulation of non-Lambertian models in the case of orthographic projection with a single light source located far from the surface. In this simplified framework, we present a unified approach to two popular models for non-Lambertian surfaces proposed by Oren-Nayar [48, 49, 47, 50] and by Phong . We will consider light sources placed in oblique directions with respect to the surface and we will use that unified formulation to develop a general numerical approximation scheme which is able to solve the corresponding nonlinear partial differential equations arising in the mathematical description of these models.
To better understand the contribution of this paper, let us start from the classical SfS problem where the goal is to reconstruct the surface from a single image. In mathematical terms, given the shading informations contained in a single two-dimensional gray level digital image , where , we look for a surface that corresponds to its shape (hence the name Shape from Shading). This problem is described in general by the image irradiance equation introduced by Bruss 
where the normalized brightness of the given grey-value image is put in relation with the function that represents the reflectance map giving the value of the light reflection on the surface as a function of its orientation (i.e., of the normal ) at each point . Depending on how we describe the function , different reflection models are determined. In the literature, the most common representation of is based on the Lambertian model (the L–model in the sequel) which takes into account only the angle between the outgoing normal to the surface and the light source , that is
denotes the standard scalar product between vectors andindicates the diffuse albedo, i.e. the diffuse reflectivity or reflecting power of a surface. It is the ratio of reflected radiation from the surface to incident radiance upon it. Its dimensionless nature is expressed as a percentage and it is measured on a scale from zero for no reflection of a perfectly black surface to 1 for perfect reflection of a white surface. The data are the grey-value image , the direction of the light source represented by the unit vector and the albedo . The light source is a unit vector, hence . In the simple case of a vertical light source, that is when the light source is in the direction of the vertical axis, this gives rise to an eikonal equation. Several questions arise, even in the simple case: is a single image sufficient to determine the surface? If not, which set of additional informations is necessary to have uniqueness? How can we compute an approximate solution? Is the approximation accurate? It is well known that for Lambertian surfaces there is no uniqueness and other informations are necessary to select a unique surface (e.g. the height at each point of local maximum for ). However, rather accurate schemes for the classical eikonal equation are now available for the approximation. Despite its simplicity, the Lambertian assumption is very strong and does not match with many real situations that is why we consider in this paper some non-Lambertian models trying to give a unified mathematical formulation for these models.
In order to set this paper into a mathematical perspective, we should mention that the pioneering work of Horn [28, 29] and his activity with his collaborators at MIT [30, 31] produced the first formulation of the Shape from Shading (SfS) problem via a partial differential equation (PDE) and a variational problem. These works have inspired many other contributions in this research area as one can see looking at the extensive list of references in the two surveys [85, 19]. Several approaches have been proposed, we can group them in two big classes (see the surveys [85, 19]): methods based on the resolution of PDEs and optimization methods based on a variational approximation. In the first group the unknown is directly the height of the surface , one can find here rather old papers based on the method of characteristics [17, 60, 29, 46, 45, 6, 38] where one typically looks for classical solutions. More recently, other contributions were developed in the framework of weak solutions in the viscosity sense starting from the seminal paper by Rouy and Tourin  and, one year later, by Lions-Rouy and Tourin  (see e.g. [33, 9, 11, 22, 10, 56, 4, 62, 36, 23, 53]). The second group contains the contribution based on minimization methods for the variational problem where the unknown are the partial derivatives of the surface, and (the so-called normal vector field. See e.g. [30, 16, 67, 26, 68, 32, 7, 82]). It is important to note that in this approach one has to couple the minimization step to compute the normal field with a local reconstruction for which is based usually on a path integration. This necessary step has also been addressed by several authors (see  and references therein). We should also mention that a continuous effort has been made by the scientific community to take into account more realistic reflectance models [3, 1, 59, 76, 77], different scenarios including perspective camera projection [44, 70, 54, 2, 75, 14] and/or multiple images of the same object [83, 84]. The images can be taken from the same point of view but with different light sources as in the photometric stereo method [81, 37, 42, 69, 43] or from different points of view but with the same light source as in stereo vision . Recent works have considered more complicated scenarios, e.g. the case when light source is not at the optical center under perspective camera projection 
. It is possible to consider in addition other supplementary issues, as the estimation of the albedo[86, 5, 65, 64] or of the direction of the light source that are usually considered known quantities for the model but in practice are hardly available for real images. The role of boundary conditions which have to be coupled with the PDE is also a hard task. Depending on what we know, the model has to be adapted leading to a calibrated or uncalibrated problem (see [84, 83, 25, 57] for more details). In this work we will assume that the albedo and the light source direction are given.
Regarding the modeling of non-Lambertian surfaces we also want to mention the important contribution of Ahmed and Farag . These authors have adopted the modeling for SfS proposed by Prados and Faugeras [55, 52] using a perspective projection where the light source is assumed to be located at the optical center of the camera instead at infinity and the light illumination is attenuated by a term ( represents here the distance between the light source and the surface). They have derived the Hamilton-Jacobi (HJ) equation corresponding to the Oren-Nayar model, developed an approximation via the Lax-Friedrichs sweeping method. They gave there an experimental evidence that the non-Lambertian model seems to resolve the classical concave/convex ambiguity in the perspective case if one includes the attenuation term . In  they extended their approach for various image conditions under orthographic and perspective projection, comparing their results for the orthographic L–model shown in  and in . Finally, we also want to mention the paper by Ragheb and Hancock  where they treat a non-Lambertian model via a variational approach, investigating the reflectance models described by Wolff and by Oren and Nayar [80, 79, 50].
In this paper we will adopt the PDE approach in the framework of weak solutions for Hamilton-Jacobi equations. As we said, we will focus our attention on a couple of non-Lambertian reflectance models: the Oren-Nayar and the Phong models [48, 49, 47, 50, 51]. Both models are considered for an orthographic projection and a single light source at infinity, so no attenuation term is considered here. We are able to write the Hamilton-Jacobi equations in the same fixed point form useful for their analysis and approximation and, using the exponential change of variable introduced by Kružkov in , we obtain natural upper bound for the solution. Moreover, we propose a semi-Lagrangian approximation scheme which can be applied to both the models, prove a convergence result for our scheme that can be applied to this class of Hamilton-Jacobi equations, hence to both non-Lambertian models. Numerical comparisons will show that our approach is more accurate also for the 3D reconstructions of non-smooth surfaces.
A similar formulation for the Lambertian SfS problem with oblique light direction has been studied in  and here is extended to non-Lambertian models. We have reported some preliminary results just for the Oren-Nayar model in .
Organization of the paper.
The paper is organized as follows. After a formulation of the general model presented in Section 2, we present the SfS models starting from the classical Lambertian model (Section 3). In Sections 4 and 5 we will give details on the construction of the nonlinear partial differential equation which corresponds respectively to the Oren-Nayar and the Phong models. Despite the differences appearing in these non-Lambertian models, we will be able to present them in a unified framework showing that the Hamilton-Jacobi equations for all the above models share a common structure. Moreover, the Hamiltonian appearing in these equations will always be convex in the gradient . Then, in Section 6, we will introduce our general approximation scheme which can be applied to solve this class of problems. In Section 7 we will apply our approximation to a series of benchmarks based on synthetic and real images. We will discuss some issue like accuracy, efficiency and the capability to obtain the maximal solution showing that the semi-Lagrangian approximation is rather effective even for real images where several parameters are unknown. Finally, in the last section we will give a summary of the contributions of this work with some final comments and future research directions.
2 Formulation of the general model
We fix a camera in a three-dimensional coordinate system (Oxyz) in such a way that Oxy coincides with the image plane and Oz with the optical axis.
Let (with ) be the unit vector that represents the direction of the light source (the vector points from the object surface to the light source);
let be the function that measures the gray-level of the input image at the point . is the datum in the model since it is measured at each pixel of the image, for example in terms of a greylevel (from to ). In order to construct a continuous model, we will assume that takes real values in the interval , defined in a compact domain called “reconstruction domain” (with open set),
, where the points with a value of are the dark point (blacks), while those with a value of correspond to a completely reflection of the light (white dots, with a maximum reflection).
We consider the following assumptions:
there is a single light source placed at infinity in the direction (the light rays are, therefore, parallel to each other);
the observer’s eye is placed at an infinite distance from the object you are looking at (i.e. there is no perspective deformation);
there are no autoreflections on the surface.
In addition to these assumptions, there are other hypothesis that depend on the different reflectance models (we will see them in the description of the individual models).
Being valid the assumption of orthographic projection, the visible part of the scene is a graph and the unit normal to the regular surface at the point corresponding to is given by:
where is the outgoing normal vector.
We assume that the surface is standing on a flat background so the height function, which is the unknown of the problem, will be non negative, . We will denote by the region inside the silhouette and we will assume (just for technical reasons) that is an open and bounded subset of (see Fig. 1).
It is well known that the SfS problem is described by the image irradiance equation (1) and depending on how we describe the function different reflection models are determined. We describe below some of them. To this end, it would be useful to introduce a representation of the brightness function where we can distinguish different terms representing the contribution of ambient, diffuse reflected and specular reflected light. We will write then
where , and are respectively the above mentioned components and , and indicate the percentages of these components such that their sum is equal to 1 (we do not consider absorption phenomena). Note that the diffuse or specular albedo is inside the definition of or , respectively. In the sequel, we will always consider normalized in . This will allow to switch on and off the different contributions depending on the model. Let us note that the ambient light term represents light everywhere in a given scene. In the whole paper we will consider it as a constant and we will neglect its contribution fixing . Moreover, for all the models presented below we will suppose uniform diffuse and/or specular albedo and we will put them equal to , that is all the points of the surface reflect completely the light that hits them. We will omit them in what follows. As we will see in the following sections, the intensity of diffusely reflected light in each direction is proportional to the cosine of the angle between surface normal and light source direction, without taking into account the point of view of the observer, but another diffuse model (the Oren–Nayar model) will consider it in addition. The amount of specular reflected light towards the viewer is proportional to , where is the angle between the ideal (mirror) reflection direction of the incoming light and the viewer direction, being a constant modelling the specularity of the material. In this way we have a more general model and, dropping the ambient and specular component, we retrieve the Lambertian reflection as a special case.
3 The Lambertian model (L–model)
For a Lambertian surface, which generates a purely diffuse model, the specular component does not exist, then in (4) we have just the diffuse component on the right side. Lambertian shading is view independent, hence the irradiance equation (1) becomes
Under these assumptions, the orthographic SfS problem consists in determining the function that satisfies the equation (5). The unit vector and the function are the only quantities known.
For Lambertian surfaces [30, 31], just considering an orthographic projection of the scene, it is possible to model the SfS problem via a nonlinear PDE of the first order which describes the relation between the surface (our unknown) and the brightness function . In fact, recalling the definition of the unit normal to a graph given in (3), we can write (5) as
where . This is an Hamilton-Jacobi type equation which does not admit in general a regular solution. It is known that the mathematical framework to describe its weak solutions is the theory of viscosity solutions as in .
The vertical light case.
If we choose , the equation (6) becomes the so-called “eikonal equation”:
The points where assumes maximum value correspond to the case in which and have the same direction: these points are usually called “singular points”.
In order to make the problem well-posed, we need to add boundary conditions to the equations (6) or (7): they can require the value of the solution (Dirichlet boundary conditions type), or the value of its normal derivative (Neumann boundary conditions), or an equation that must be satisfied on the boundary (the so-called boundary conditions “state constraint”). In this paper, we consider Dirichlet boundary conditions equal to zero assuming a surface on a flat background
but a second possibility of the same type occurs when it is known the value of on the boundary, which leads to the more general condition
Unfortunately, adding a boundary condition to the PDE that describes the SfS model is not enough to obtain a unique solution because of the concave/convex ambiguity. In fact, the Dirichlet problem (6)-(10) can have several weak solutions in the viscosity sense and also several classical solutions due to this ambiguity (see ). As an example, all the surfaces represented in Fig. 2 are viscosity solutions of the same problem (7)-(9) which is a particular case of (6)-(10) (in fact the equation is with homogenous Dirichlet boundary condition). The solution represented in Fig. 2-a is the maximal solution and is smooth. All the non-smooth a.e. solutions, which can be obtained by a reflection with respect to a horizontal axis, are still admissible weak solutions (see Fig. 2-b). In this example, the lack of uniqueness of the viscosity solution is due to the existence of a singular point where the right hand side of (7) vanishes. An additional effort is then needed to define which is the preferable solution since the lack of uniqueness is also a big drawback when trying to compute a numerical solution. In order to circumvent these difficulties, the problem is usually solved by adding some information such as height at each singular point .
For analytical and numerical reasons it is useful to introduce the exponential Kružkov change of variable  . In fact, setting the problem in the new variable we will have values in instead of as the original variable so an upper bound will be easy to find. Note that is a free positive parameter which does not have a specific physical meaning in the SfS problem. However, it can play an important role also in our convergence proof as we will see later (see the remark following the end of Theorem 6.1). Assuming that the surface is standing on a flat background and following , we can write (6) and (9) in a fixed point form in the new variable . To this end let us define and as
and let denote the unit ball in . We obtain
It is important to note for the sequel that the structure of the above first order Hamilton-Jacobi equation is similar to that related to the dynamic programming approach in control theory where is a vector field describing the dynamics of the system and is a running cost. In that framework the meaning of is that of a value function which allows to characterize the optimal trajectories (here they play the role of characteristic curves). The interested reader can find more details on this interpretation in .
4 The Oren-Nayar model (ON–model)
The diffuse reflectance ON–model [48, 49, 47, 50] is an extension of the previous L–model which
explicitly allows to handle rough surfaces.
The idea of this model is to represent a rough surface as an aggregation of V-shaped cavities, each with Lambertian reflectance properties (see Fig. 3).
In  and, with more details, in , Oren and Nayar derive a reflectance model for several type of surfaces with different slope-area distributions. In this paper we will refer to the model called by the authors the “Qualitative Model”, a simpler version obtained by ignoring interreflections (see Section of  for more details).
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
Note that and are two non negative constants depending on the statistics of the cavities via the roughness parameter . We set , interpreting as the slope of the cavities. In this model (see Fig. 4), represents the angle between the unit normal to the surface and the light source direction , stands for the angle between and the observer direction , is the angle between the projection of the light source direction and the axis onto the -plane, denotes the angle between the projection of the observer direction and the axis onto the -plane and the two variables and are given by
Since the vectors and are fixed and given, their projection on the incident plane is obtained considering their first two components over three (see Eq. (21)). In this way, the quantity is computed only once for a whole image.
We define (see Fig. 4):
For smooth surfaces, we have and in this case the ON–model reduces to the L–model. In the particular case , or, more precisely, when (e.g. the case when the unit vectors and are perpendicular we get ) the equation (14) simplifies and reduces to a L–model scaled by the coefficient . This means that the model is more general and flexible than the L–model. This happens when only one of the two unit vectors is zero or, more in general, when the dot product between the normalized projections onto the -plane of and is equal to zero.
Also for this diffuse model we neglect the ambient component, setting . As a consequence, in the general equation (4) the total light intensity is equal to the diffuse component (described by the equation (14)). This is why we write instead of in what follows.
To deal with this equation one has to compute the and operators which appear in (14) and (18).
Hence, we must consider several cases described in detail in what follows. For each case we will derive a partial differential equation that is always a first order nonlinear HJ equation:
Case 1: and
The brightness equation (14) becomes
where and .
Case 2: and
In this case the brightness equation (14) becomes
Case 3: and
In this case we have the implication . The brightness equation (14) simplifies in
and the HJ equation associated to it becomes consequentially
that is equal to the L–model scaled by the coefficient .
Case 4: and
This is a particular case when the position of the light source coincides with the observer direction but there are not on the vertical axis. This choice implies , then defining , the equation (14) simplifies to
and we arrive to the following HJ equation
Note that this four cases are exactly the same cases reported and analyzed in . This is not surprising since the reflectance model used there is always the same one proposed by Oren and Nayar. However, here we get different HJ equations since we consider an orthographic camera projection and cartesian coordinates whereas in  the HJ equations are derived in spherical coordinates under a generalized perspective camera projection. Another major difference is that in that paper the light source is close to the camera but is not located at the optical center of the camera.
The vertical light case.
If , independently of the position of , the analysis is more simple. In fact, the first three cases considered above reduce to a single case corresponding to the following simplified PDE for the brightness equation (14)
In this way we can put it in the following eikonal type equation, analogous to the Lambertian eikonal equation (7):
using the equivalence
Let us define the vector field for the ON-model
Then, introducing the exponential Kružkov change of variable as already done for the L–model, we can finally write the fixed point problem in the new variable obtaining the
Note that the simple homogeneous Dirichlet boundary condition is due to the flat background behind the object but a condition like can also be considered if necessary. Moreover, the structure is similar to the previous Lambertian model although the definition of the vector field and of the cost are different.
In the particular case when , the equation (14) simply reduces to
and, as a consequence, the Dirichlet problem in the variable is equal to (42) with .
5 The Phong model (PH–model)
The PH–model is an empirical model that was developed by Phong  in 1975. This model introduces a specular component to the brightness function , representing the diffuse component in (4) as the Lambertian reflectance model.
A simple specular model is obtained putting the incidence angle equal to the reflection one and , and belong to the same plane.
This model describes the specular light component as a power of the cosine of the angle between the unit vectors and (it is the vector representing the reflection of the light on the surface). Hence, the brightness equation for the PH–model is
where the parameter is used to express the specular reflection properties of a material and and indicate the percentages of diffuse and specular components, respectively. Note that the contribution of the specular part decreases as the value of increases.
Starting to see in details the PH–model in the case of oblique light source and oblique observer .
Assuming that is the bisector of the angle between and (see Fig. 5), we obtain
From the parallelogram law, taking into account that and are unit vectors, we can write , then we can derive the unit vector as follow:
With this definition of the unit vector we have
Then, setting , that represents the worst case, equation (44) becomes
to which we add a Dirichlet boundary condition equal to zero, assuming that the surface is standind on a flat blackground.
Note that the cosine in the specular term is usually replaced by zero if (and in that case we get back to the L–model).
As we have done for the previous models, we write the surface as , for , , and , so (5) will be written as
Dividing by , defining as in (36) and , we get
By the equivalence we obtain
Let us define the vector field
Let us also define
Again, using the exponential Kružkov change of variable as done for the previous models, we can finally write the nonlinear fixed point problem
A. Oblique light source and vertical position of the observer.
In the case of oblique light source and vertical observer , the dot product becomes
The fixed point problem in will be equal to (57) with the following choices
B. Vertical light source and oblique position of the observer.
When the definition of the vector reported in (5) becomes
and, as a consequence, the dot product with general is
Hence, the fixed point problem in is equal to (57) with
C. Vertical light source and vertical position of the observer.
If we choose the equation (5) simplifies in
Working on this equation one can put it in the following eikonal type form, which is analogous to the Lambertian eikonal equation (7):
Remark on the control interpretation. The above analysis has shown that all the cases corresponding to the models proposed by Oren-Nayar and by Phong lead to a stationary Hamilton-Jacobi equation of the same form, namely
where the vector field and the cost can vary according to the model and to the case. This gives to these models a control theoretical interpretation which can be seen as a generalization of the control interpretation for the original Lambertian model (which was related to the minimum time problem). In this framework, is the value function of a rescaled (by the Kružkov change of variable) control problem in which one wants to drive the controlled system governed by
( here is the control function taking values in ) from the initial position to the target (the silhouette of the object) minimizing the cost associated to the trajectory. The running cost associated to the position and the choice of the control will be given by . More informations on this interpretation, which is not crucial to understand the application to the SfS problem presented in this paper, can be found in .
6 Semi-Lagrangian Approximation
Now, let us state a general convergence theorem suitable for the class of differential operators appearing in the models described in the previous sections. As we noticed, the unified approach presented in this paper has the big advantage to give a unique formulation for the three models in the form of a fixed point problem
where indicates the model, i.e. .
We will see that the discrete operators of the ON–model and the PH–model described in the previous sections satisfy the properties listed here. In order to obtain the fully discrete approximation we will adopt the semi-Lagrangian approach described in the book by Falcone and Ferretti . The reader can also refer to  for a similar approach to various image processing problems (including nonlinear diffusions models and segmentation).
Let so that will be the vector solution giving the approximation of the height at every node of the grid. Note that in one dimension the index is an integer number, in two dimensions denotes a multi-index, . We consider a semi-Lagrangian scheme written in a fixed point form, so we will write the fully discrete scheme as
Denoting by the global number of nodes in the grid, the operator corresponding to the oblique light source is that is defined componentwise by
represents an interpolation operator based on the values at the grid nodes and
Since is approximated via by interpolation on (which is defined on the grid ), it is important to use a monotone interpolation in order to preserve the properties of the continuous operator in the discretization. To this end, the typical choice is to apply a piecewise linear (or bilinear) interpolation operator which allows to define a function defined for every (and not only on the nodes)
A simple explanation for (76)-(77) is that the coefficients represent the local coordinates of the point with respect to the grid nodes (see  for more details and other choices of interpolation operators). Clearly, in (71) we will apply the interpolation operator to the point and we will denote by the function defined by .
Comparing (71) with its analogue for the vertical light case we can immediately note that the former has the additional term which requires analysis.
Let the -th component of the operator defined as in (71). Then, the following properties hold true:
Then implies .
is a monotone operator, i.e., implies .
is a contraction mapping in if
To prove that implies we just note that
Let ; then
This implies that if