A Unified Shape-From-Shading Approach for 3D Surface Reconstruction Using Fast Eikonal Solvers

. Object shape reconstruction from images has been an active topic in computer vision. Shape-from-shading (SFS) is an important approach for inferring 3D surface from a single shading image. In this paper, we present a uniﬁed SFS approach for surfaces of various reﬂectance properties using fast eikonal solvers. The whole approach consists of three main components: a uniﬁed SFS model, a uniﬁed eikonal-type partial diﬀerential image irradiance (PDII) equation, and fast eikonal solvers for the PDII equation. The ﬁrst component is designed to address diﬀerent reﬂectance properties including diﬀuse, specular, and hybrid reﬂections in the imaging process of the camera. The second component is meant to derive the PDII equation under an orthographic camera projection and a single distant point light source whose direction is the same as the camera. Finally, the last component is targeted at solving the resultant PDII equation by using fast eikonal solvers. It comprises two Godunov-based schemes with fast sweeping method that can handle the eikonal-type PDII equation. Experiments on several synthetic and real images demonstrate that each type of the surfaces can be eﬀectively reconstructed with more accurate results and less CPU running time.


Introduction
In the field of computer vision, object shape reconstruction from images has been an active topic.
ere are several techniques, such as stereo vision, structured light, fringe projection profilometry, and shape-from-X (X � shading, photometric stereo, texture, focus/defocus, motion, etc.). Shape-from-shading (SFS) is an important approach for inferring 3D surface from a single shading image and because of its simplicity of equipment, it is widely used in face reconstruction [1,2], 3D reconstruction of medical images [3,4], lunar surface reconstruction [5,6], and so on. It was initiated by Horn [5] who firstly formulated a first-order partial differential image irradiance (PDII) equation describing the relations between the 3D shape of a surface and its corresponding 2D variation of intensities. us one can determine 3D surface by starting with the PDII equation.
Since Horn's original work, a great number of different SFS approaches have come out (for surveys, refer to Zhang et al. [7], and Durou et al. [8]).
ere are mainly two steps when utilizing an SFS approach. e first step is meant to model the image formation process of the camera which is determined by the reflectance property of the surface, the light source, and the camera projection and to derive the PDII equation under certain assumption [9]. e second step is targeted at designing a numerical scheme to solve the resultant PDII equation. Most of the SFS approaches concentrate on how to design an effective numerical scheme assuming that the surface obeys a simple Lambertian reflection.
ese approaches are generally divided into two classes: partial differential equation-(PDE-) based methods and optimisation methods [10,11]. e characteristics-based approach [5] and the viscosity solution-based approaches [1,3,4,[9][10][11][12][13][14][15][16] can be categorized into the first class. We should mention the pioneering viscosity solution-based approach of Rouy and Tourin [12], who first described the PDII equation under Lambertian reflectance model as a Hamilton-Jacobi-Bellman PDE and got a nonclassical solution based on viscosity solution theory. Kimmel and Sethian [13] transformed the PDII equation under the vertical light into an eikonal-type PDE and used the firstorder fast marching method [17] to solve its viscosity solution. For the oblique light case, Governi et al. [14] reconstructed the initial surface by directly using the fast marching method [17]. ey rotated the normal map obtained from the surface around the oblique light and then computed the "new" image as the dot product between the normal map and vertical light. e final surface could be reconstructed by applying the fast marching method to the "new" image again. e works of [5,7,8,[10][11][12][13][14] are thinking of an orthographic projection for the camera. As for perspective camera projection, Prados and Faugeras [1] related the PDII equation to a Hamiltonian based on the work [12] and got its viscosity solution with optimal control theory. Breuss et al. [15,16] analytically and numerically studied the perspective PDII equation formulated by Prados and Faugeras [1] and the associated Hamilton-Jacobi PDE. At the same time, they proved the convergence of the finite-difference and the semi-Lagrangian schemes for the resultant PDE. e second class means the minimisation methods for the SFS problem [7]. Ikeuchi and Horn [18] formulated SFS as a minimisation problem of the difference between observed intensities and the expected intensities that are given through the PDII equation from the expected surface normal, on which the smoothness constraint was used. Tankus et al. [19] first derived a perspective PDII equation and obtained an approximate solution under the assumption that the surface is locally paraboloidal. e 3D shape was reconstructed by minimising a quadratic cost functional. More recently, Santo et al. [20] revisited the numerical SFS approach of Ikeuchi and Horn [18] and described corresponding solution that was built upon different convex relaxation strategies. It is worth mentioning that Quéau et al. [21] combined the advantages of optimisation methods and PDE-based methods and built a generic variational solution that is suitable for SFS under natural illumination and can handle a variety of scenarios for various lighting and camera projection.
While most of the SFS approaches assume the Lambertian reflection, there are a few researchers who are interested in non-Lambertian SFS since the Lambertian reflectance model has been proved to be inaccurate, especially for rough diffuse surfaces [22]. Ragheb and Hancock [23] proposed a non-Lambertian SFS with the Oren-Nayar reflectance model and gave two solutions: the lookup table and the analytic solution. Ahmed and Farag [24,25] presented several non-Lambertian SFS approaches including Ward SFS and Oren-Nayar SFS and approximated the PDII equations by using the Lax-Friedrichs sweeping scheme [26]. Since the actual convergence to the correct solution is very slow in [25], Vogel and Cristiani [27] applied the Upwind scheme to get a more efficient solution with less convergence time. Tozza and Falcone [10,28] addressed a general framework for several non-Lambertian SFS problems including Oren-Nayar SFS and Phong SFS, solved by a semi-Lagrangian scheme, and obtained convergence results. However, their framework can only handle a special case where the specular reflection parameter n in the Phong reflectance model [29] equals 1; that is, it represents the worst case. By extending the work of Galliani et al. [30], Ju et al. [4] adopted spherical parameterisation of the surface into the Oren-Nayar PDII equations and thus could compute them at any position of the point light source. However, the fast-marching scheme depicted in Cartesian coordinates needs to be converted to spherical coordinates during the process.
In this paper, motivated by the work of Camilli and Tozza [31] and based on our previous work [11], we first present a unified SFS model for surfaces of different reflectance properties including diffuse, specular, and hybrid reflections in the image formation process. Although our work falls in the situation where the camera performs an orthographic projection and the direction of the single distant point light source is the same as the camera, these reflections lead to more complex nonlinear PDII equations. However, all the PDII equations corresponding to the reflections considered here (Oren-Nayar model and unified model) have a similar structure, so we can look for weak solutions to this class in the viscosity solution sense. Another contribution of our work is that we convert the PDII equation into an eikonal-type PDE through solving a highorder equation by using the Newton-Raphson method, after which we try to obtain the viscosity solution of the eikonaltype PDE by using fast eikonal solvers which are composed of the first-and high-order Godunov-based schemes accelerated by the fast sweeping method.
A similar formulation for the SFS problem of the Oren-Nayar model has been reported in our previous work [11]. As we said, in this paper we will focus our attention on the unified reflectance model (including Lambertian model, Oren-Nayar model, and Blinn-Phong model) and formulate the unified high-order PDII equation under the vertical light. Using the Newton-Raphson method for the resultant PDII equation, we will obtain the eikonal-type PDE that can be solved via fast eikonal solvers presented preliminarily in our work [11].

A Unified SFS Model in the Imaging Process
In this section, a very brief description for the Lambertian, Oren-Nayar and Blinn-Phong reflectance model is given in order to setup a unified imaging model.

Lambertian Reflectance Model.
Generally, the Lambertian reflectance is a classical assumption in most of the SFS approaches [1, 3, 5, 8, 12-16, 18-21, 30] for approximating the reflectance property of the diffuse surface. In this case, the surface reflected radiance is addressed as [3] where I 0 is the intensity of point light source, ρ d is the diffuse albedo which controls the proportion of incident light that is reflected diffusely, and θ i is the angle between the surface unit normal n and the incident light direction L illustrated in Figure 1.

Oren-Nayar Reflectance Model.
In order to get rid of the inaccuracy resulting from the assumption of the Lambertian reflectance model for diffuse reflection, Oren and Nayar [22] proposed a comprehensive reflectance model for rough diffuse surfaces. By assuming that the surface is composed of V-shaped cavities which are symmetric and have two planar facets and that each facet obeys a simple Lambertian reflection, for a Gaussian distribution of the facet normals, they got a simplified expression for the reflected radiance: e parameter σ is applied as a measure of the surface roughness, and it denotes the standard deviation of the Gaussian distribution.
For smooth surfaces, we have σ � 0 and obviously the Oren-Nayar reflectance model degenerate to the Lambertian model in this situation.

Blinn-Phong Reflectance Model.
It is worth mentioning that Phong [29] developed a hybrid reflectance model by introducing a specular component to the surface reflected radiance (1). He described the specular component as a power of the cosine of the angle between the reflected light direction R and the camera direction V. Hence, the hybrid reflected radiance can be derived in general as where w d and w s are the weighting factors of diffuse and specular components, respectively, and w d + w s ≤ 1. ρ s is the specular albedo that determines the proportion of incident light that is reflected specularly. e parameter n is used to express the specular reflection property of a surface and can be used as a measure of the surface shininess. Obviously, the contribution of the specular component decreases when the value of parameter n increases. Note that it is not convenient to compute the specular reflected radiance in terms of (R · V). e Blinn-Phong reflectance model, proposed by Blinn [32], is a modification of the Phong model for computation convenience. (3), the hybrid reflected radiance based on the Blinn-Phong model can be formulated as

A Unified Reflectance Model.
As mentioned before, the Lambertian reflectance model has been proved to be inaccurate, especially for rough diffuse surfaces. us, we can combine diffuse and specular components of a surface through a linear combination of Oren-Nayar model and the specular part of Blinn-Phong model; that is, we substitute surface reflected radiance (2) into the diffuse part of Blinn-Phong model: Obviously, surface reflected radiance (5) is a unified reflectance model including the Lambertian, the Oren-Nayar, and the Blinn-Phong model. For w s � 0, it reduces to the Oren-Nayar model. For σ � 0, it reduces to the Blinn-Phong model. Specially, if w s � 0 and σ � 0, it degenerates to the Lambertian model. e following relationship between the surface reflected radiance L r and the image irradiance E i of the camera is well known [9]: where D is the entrance pupil diameter of the camera lens whose focal length is f. χ is the angle between the line of sight to an image point of a corresponding surface point and the optical axis of the camera. Even for uniform illumination, the term cos 4 χ implies nonuniform image irradiance. e actual imaging lens of the camera, however, is generally designed to correct it. us, one can consider E i to be proportional to L r : Figure 1: Reflection geometry of a local surface point. n is the unit normal of the surface point P; (θ r , ϕ r ) and (θ i , ϕ i ) are the camera direction V and incident light direction L, respectively; R is reflected light direction; h is the unit angular bisector of V and L; that is, h � (V + L)/‖V + L‖; δ is the angle between n and h.
International Journal of Optics 3 Substituting equation (5) into (7), and if we set ρ d � ρ s to a constant ρ as done in [31] and denote I � πE i /ηI 0 ρ as done in most of SFS approaches, the image irradiance equation (7) will be rewritten as

A Unified Eikonal-Type PDII Equation
In this section, we will formulate the image irradiance equation under the situation where the camera performs an orthographic projection and the direction of the single distant point light source coincides with the camera.

Nonlinear PDII Equation for the Unified Model.
With the basis that the optical axis of the camera is the z−axis and the image plane of the camera is the x − y plane, the SFS approach can be described as inferring a 3D surface, z(x, y). Since our work falls in an orthographic camera projection, the first partial derivatives of the surface z(x, y) with respect to x and y, respectively, are So the unit normal n at a 3D surface point P(x, y, z(x, y)) can be expressed as Considering that the direction of the distant point light source L is the same as the camera direction V illustrated in Figure 1, we have θ i � θ r , ϕ i � ϕ r , α � θ i � β, and h � L/‖L‖, n · h� cos θ i . Consequently, image irradiance equation (8) will be reduced to Defining that the direction vectors of L and V both are [0, 0, −1]; that is, they are parallel to the optical axis of the camera lens, and because θ i is the angle between n and L, we have Substituting equation (12) into (11), the image irradiance equation (11) can be rewritten as Obviously, the image irradiance equation (13) is a more complex nonlinear PDE and is difficult to solve z(x, y).
Applying the change of variable T � 1/ ������������� 1 + ‖∇z(x, y)‖ 2 , the PDII equation (13) can be considered as calculating a zero of the function F(T), given by 3.2. Eikonal-Type PDE for the Oren-Nayar Model. Especially, for w s � 0, that is, for the PDII equation of the Oren-Nayar model, F(T) � 0 will lead to a quadratic equation: Calculating equation (15) and satisfying 0 < T ≤ 1, we can obtain Hence, SFS problem (13) can be rewritten as an eikonaltype PDE: where Ω is a given image domain and Γ(x, y) is a boundary condition. Similar work has been studied in our previous work [11] and here will be extended to the unified model. International Journal of Optics

Eikonal-
If the surface roughness 0 ≤ σ ≤ 0.6220, then A ≥ 2B ≥ 0. At the same time, 0 < T ≤ 1, so F ′ (T) > 0 and function F(T) is monotonous. Simultaneously, we have Hence, function (14) always has a unique zero. Starting with the value T 0 � 0, the iterative equation of the Newton-Raphson method is applied to calculate a new value for T k as follows: After several numbers of iterations, an accurate zero of function (14) is obtained. Similar to the structure of the Oren-Nayar model, we can get an eikonal-type PDE for the unified model:

Fast Eikonal Solvers for the Eikonal-Type PDE
In this section, we will use the fast eikonal solvers which are composed of the first-order Godunov-based scheme [12,33] and high-order Godunov-based scheme [11,34] accelerated by the fast sweeping method [33,35] to look for the weak solutions of the resultant eikonal-type PDE (21) in the viscosity solution sense.

First-Order Godunov-Based Scheme.
We use (x i , y j ) � (i × w, j × w) to denote a grid point in the image domain Ω, w to denote the grid size, M × N to denote the image size, and z i,j � z(x i , y j ) to denote the numerical solution at the 3D surface z(x, y). e first-order Godunov-based scheme [12,33] can be employed to discretize resultant eikonal-type PDE (21): where us, the viscosity solution of eikonal-type PDE (21) can be obtained using the first-order Godunov-based scheme:

High-Order Godunov-Based Scheme.
In order to obtain a higher-order accuracy viscosity solution, the high-order Godunov-based scheme [34] can be employed to discretize resultant eikonal-type PDE (21): where p i,j and q i,j need to be approximated with higherorder accuracy. According to [34], third-order weighted essentially nonoscillatory scheme [36] is able to be chosen as p i,j and q i,j approximations: with International Journal of Optics where ε is a very small number that keeps the denominator from getting too close to zero. Similarly, q + i,j and q − i,j can be defined. Now the viscosity solution of eikonal-type PDE (21) can be obtained using the high-order Godunov-based scheme:

Fast Sweeping Method for Godunov-Based Schemes.
In order to speed up the convergence numerical schemes, we take the philosophy of fast sweeping method [33,35] to the first-order or high-order Godunov-based schemes in the following. When the derivatives p + i,j , q + i,j and p − i,j , q − i,j are calculated, the newest available values for z are employed. Meanwhile, the iterations do not sweep in only one direction but in four alternating directions repeatedly: (1) from upper left to lower right, that is, i � 1: I, j � 1:J; (2) from lower left to upper right, that is, i � I: 1, j � 1: J; (3) from lower right to upper left, that is, i � I: 1 , j � J: 1; (4) from upper right to lower left, that is, i � 1: I, j � J: 1. As can be easily seen, various values z i±1,j , z i±2,j and z i,j±1 , z i,j±2 are to be taken according to the current sweeping direction.
We summarize the fast eikonal solvers for the resultant eikonal-type PDE (21) as follows: Step 1 (Initialization): according to the boundary condition z(x, y) � Γ(x, y), (x, y) ∈zΩ, assign exact values at the grid points on the boundary zΩ, whose values are fixed during iterations. At all other grid points, for firstorder Godunov-based scheme, big positive values are used as the initial guess, which are larger than the maximum of the true solutions and will be updated in the process of iterations. Especially for high-order Godunovbased scheme, the solution of the first-order Godunovbased scheme is considered as the initial guess.
Step 2 (Alternating Sweepings): we compute z new i,j according to the update formulation (24)  Step 3 (Convergence): if ‖z new − z old ‖ L 1 ≤ μ, where μ is a given threshold value, the schemes converge and stop; otherwise, return to Step 2. In this paper, we use μ � 10 − 5 .

Experimental Results
Several experiments on synthetic and real images with different reflectance properties have been carried out in order to assess the effectiveness of the presented unified SFS approach. We compare our presented approach with the Ahmed and Farag's approach [24,25] using Lax-Friedrichs sweeping scheme for the same reflectance property. We implement all the approaches in Matlab. All the experiments are conducted on a PC with a Xeon E5-1650 processor and 16 GB of DDR3 memory.

Experimental
Results on Synthetic Images. We use two synthetic surfaces including a ball and a vase, which have been benchmark test surfaces and are determined by equations (30) and (31), respectively: where R � 75 is the radius of the ball and the generated image size M × N is 256 × 256; that is, (x, y) ∈ [−127, 128] × [−127, 128]: where g(x) � 0.15 − 0.025(6x − 1)(2x − 1) 2 (3x + 2) 2 (2x + 1) and original range of (x, y) values is In order to assess the effectiveness of the presented unified SFS approach for the surfaces of various reflectance properties, four different parameter sets of σ, w d , w s , and n are used to generate the shading images. Table 1 shows the parameter values. Especially for set (1) and set (2), σ � 0 means that the unified model reduces to the Blinn-Phong model with different diffuse and specular components, and for set (3), it means that the unified model reduces to the Oren-Nayar diffuse model. e experimental results for the synthetic ball images are illustrated in Figure 3. Figures 3(a)-3(d) show the shading images generated by the four parameter sets shown in Table 1 Figure 4 illustrates the corresponding experimental results for the synthetic vase images.
As can be roughly seen from Figures 3 and 4, the fast eikonal solvers and the Lax-Friedrichs sweeping scheme can basically get satisfactory reconstructed results for the four different parameter sets of the unified reflectance model. Furthermore, we can easily see that the first-and high-order Godunov-based schemes illustrate similar results, and both schemes can give much better reconstructed results with smaller differences between reconstructed surfaces and ground truths than the Lax-Friedrichs sweeping scheme, especially for more specular components such as Figures 3(b), 3(d), 4(b), and 4(d).
e effectiveness of our presented unified approach is further described by comparisons between the fast eikonal solvers and the Lax-Friedrichs sweeping scheme with the mean absolute (MA) error, the root mean square (RMS) error, and the CPU running time. Tables 2 and 3 list the quantitative comparisons of the three schemes for the synthetic ball and vase images. It can be seen obviously that the first-order Godunov-based scheme shows much more superiority in CPU running time in all the images that we carried out since it converges after about 2 iterations. At the same time, we can see that the high-order Godunov-based scheme exhibits the minimal reconstructed error in both the MA and RMS errors because the third-order weighted essentially nonoscillatory scheme is adopted in the approximation process. e Lax-Friedrichs sweeping scheme shows a worse performance; maybe, it is difficult to look for a perfect estimate for the artificial viscosity term.

Experimental Results on Real Images.
In order to demonstrate the performance of our presented approach for real surface, we test it on two real images and also compare the reconstructed results with the Lax-Friedrichs sweeping scheme. e first image is a vase applied in [7], which is illustrated in Figure 5(a) and is mostly diffuse. e second image is a plastic bottle, which is illustrated in Figure 6 (e), we can draw the similar conclusions as for more specular surface. As shown previously, from Figures 5(e) and 6(e), we can see that the Lax-Friedrichs sweeping scheme also exhibits a slightly worse performance since it is hard to find a perfect estimate for the artificial viscosity term. It is well worth noting that the firstorder Godunov-based scheme is the fastest and the reconstructed surface using the high-order Godunov-based scheme looks like the best result. International Journal of Optics

Conclusions
In this paper, we have reported a unified SFS approach for surfaces of various reflectance properties including diffuse, specular, and hybrid reflections using fast eikonal solvers. A unified reflectance model that is a linear combination of the Oren-Nayar model and the specular part of the Blinn-Phong model is presented. We have derived the unified image irradiance equation under this unified model with an orthographic camera projection and a single distant point light source whose direction is the same as the camera. We have also converted the PDII equation into an eikonal-type PDE through solving a high-order equation by using the Newton-Raphson method. Fast eikonal solvers which are comprised of the first-and high-order Godunov-based schemes accelerated by the fast sweeping method are employed to solve the viscosity solution of the resultant eikonal-type PDE. Finally, the experiments are conducted on both synthetic and real images and the results verify that our presented approach can provide satisfactory 3D surface reconstruction with a higher accuracy in less CPU running time.
Frankly speaking, the presented unified SFS approach can only handle the special case which assumes an orthographic camera projection and a single distant point light source whose direction is parallel to the optical axis of the camera lens. In future work, we will adopt the idea of using the Newton-Raphson method to solve the high-order PDII equations derived from the SFS problem with a more complex reflectance model and will relax the two assumptions by employing a nearby point light source and a perspective camera projection. e attenuation term of the light illumination will be also considered to eliminate the convexconcave ambiguity which can make the SFS problem illposed.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.