Abstract
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 unified SFS approach for surfaces of various reflectance properties using fast eikonal solvers. The whole approach consists of three main components: a unified SFS model, a unified eikonal-type partial differential image irradiance (PDII) equation, and fast eikonal solvers for the PDII equation. The first component is designed to address different reflectance properties including diffuse, specular, and hybrid reflections 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 effectively reconstructed with more accurate results and less CPU running time.
1. Introduction
In the field of computer vision, object shape reconstruction from images has been an active topic. There 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. Thus 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]). There are mainly two steps when utilizing an SFS approach. The 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]. The 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. These approaches are generally divided into two classes: partial differential equation- (PDE-) based methods and optimisation methods [10, 11]. The characteristics-based approach [5] and the viscosity solution-based approaches [1, 3, 4, 9–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 first-order 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]. They 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. The final surface could be reconstructed by applying the fast marching method to the “new” image again. The works of [5, 7, 8, 10–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. The 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. The 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 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 high-order equation by using the Newton-Raphson method, after which we try to obtain the viscosity solution of the eikonal-type 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].
2. 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.
2.1. 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 is the intensity of point light source, is the diffuse albedo which controls the proportion of incident light that is reflected diffusely, and is the angle between the surface unit normal and the incident light direction illustrated in Figure 1.
2.2. 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:where , ; is the incident light direction ; is the camera direction ; , . The 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 and obviously the Oren-Nayar reflectance model degenerate to the Lambertian model in this situation.
2.3. 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 and the camera direction . Hence, the hybrid reflected radiance can be derived in general aswhere and are the weighting factors of diffuse and specular components, respectively, and . is the specular albedo that determines the proportion of incident light that is reflected specularly. The parameter 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 increases.
Note that it is not convenient to compute the specular reflected radiance in terms of . The Blinn–Phong reflectance model, proposed by Blinn [32], is a modification of the Phong model for computation convenience. Substituting into in formula (3), the hybrid reflected radiance based on the Blinn–Phong model can be formulated as
2.4. A Unified Reflectance Model
As mentioned before, the Lambertian reflectance model has been proved to be inaccurate, especially for rough diffuse surfaces. Thus, 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 , it reduces to the Oren–Nayar model. For , it reduces to the Blinn–Phong model. Specially, if and , it degenerates to the Lambertian model.
The following relationship between the surface reflected radiance and the image irradiance of the camera is well known [9]:where is the entrance pupil diameter of the camera lens whose focal length is . 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 implies nonuniform image irradiance. The actual imaging lens of the camera, however, is generally designed to correct it. Thus, one can consider to be proportional to :
Substituting equation (5) into (7), and if we set to a constant as done in [31] and denote as done in most of SFS approaches, the image irradiance equation (7) will be rewritten as
3. 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.
3.1. Nonlinear PDII Equation for the Unified Model
With the basis that the optical axis of the camera is the axis and the image plane of the camera is the plane, the SFS approach can be described as inferring a 3D surface, . Since our work falls in an orthographic camera projection, the first partial derivatives of the surface with respect to and , respectively, are
So the unit normal at a 3D surface point can be expressed as
Considering that the direction of the distant point light source is the same as the camera direction illustrated in Figure 1, we have , and , . Consequently, image irradiance equation (8) will be reduced to
Defining that the direction vectors of and both are ; that is, they are parallel to the optical axis of the camera lens, and because is the angle between and , 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 . Applying the change of variable , the PDII equation (13) can be considered as calculating a zero of the function , given by
3.2. Eikonal-Type PDE for the Oren–Nayar Model
Especially, for , that is, for the PDII equation of the Oren–Nayar model, will lead to a quadratic equation:
Calculating equation (15) and satisfying , we can obtain
Hence, SFS problem (13) can be rewritten as an eikonal-type PDE:where is a given image domain and is a boundary condition. Similar work has been studied in our previous work [11] and here will be extended to the unified model.
3.3. Eikonal-Type PDE for the Unified Model
For , is a high-order function of when and it is difficult to calculate the zero values. We can use the Newton–Raphson method to solve it. The derivative of the function is
If the surface roughness , then . At the same time, , so and function is monotonous. Simultaneously, we have
Hence, function (14) always has a unique zero. Starting with the value , the iterative equation of the Newton-Raphson method is applied to calculate a new value for 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:
4. 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.
4.1. First-Order Godunov-Based Scheme
We use to denote a grid point in the image domain , to denote the grid size, to denote the image size, and to denote the numerical solution at the 3D surface . The first-order Godunov-based scheme [12, 33] can be employed to discretize resultant eikonal-type PDE (21):where , , , and
Thus, the viscosity solution of eikonal-type PDE (21) can be obtained using the first-order Godunov-based scheme:
4.2. 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):withwhere and need to be approximated with higher-order accuracy. According to [34], third-order weighted essentially nonoscillatory scheme [36] is able to be chosen as and approximations:withwhere is a very small number that keeps the denominator from getting too close to zero. Similarly, and can be defined. Now the viscosity solution of eikonal-type PDE (21) can be obtained using the high-order Godunov-based scheme:
4.3. 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 , and , are calculated, the newest available values for 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, :J; (2) from lower left to upper right, that is, ; (3) from lower right to upper left, that is, ; (4) from upper right to lower left, that is, . As can be easily seen, various values , and , 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 , assign exact values at the grid points on the boundary , whose values are fixed during iterations. At all other grid points, for first-order 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 Godunov-based scheme, the solution of the first-order Godunov-based scheme is considered as the initial guess. Step 2 (Alternating Sweepings): we compute according to the update formulation (24) or (29) by Gauss-Seidel iterations with four alternating direction sweepings:(1);(2);(3);(4). Step 3 (Convergence): if , where is a given threshold value, the schemes converge and stop; otherwise, return to Step 2. In this paper, we use .
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.
5.1. 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 is the radius of the ball and the generated image size is ; that is, :where and original range of values is . To obtain the same image size as the ball, we map range to and scale simultaneously. Their ground truths are shown in Figure 2.
(a)
(b)
In order to assess the effectiveness of the presented unified SFS approach for the surfaces of various reflectance properties, four different parameter sets of , and are used to generate the shading images. Table 1 shows the parameter values. Especially for set (1) and set (2), 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.
The 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, respectively. Figures 3(e)–3(h) show the reconstructed surfaces of Figures 3(a)–3(d) using the first-order Godunov-based scheme, while Figures 3(i)–3(l) show the reconstructed surfaces using the high-order Godunov-based scheme. Finally, Figures 3(m)–3(p) show the reconstructed surfaces of Figures 3(a)–3(d) using the Lax-Friedrichs sweeping scheme. Figure 4 illustrates the corresponding experimental results for the synthetic vase images.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
(p)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
(p)
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).
The 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. The Lax-Friedrichs sweeping scheme shows a worse performance; maybe, it is difficult to look for a perfect estimate for the artificial viscosity term.
5.2. 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. The first image is a vase applied in [7], which is illustrated in Figure 5(a) and is mostly diffuse. The second image is a plastic bottle, which is illustrated in Figure 6(a) and contains more specular components. Figures 5(b) and 6(b) show the masks of Figures 5(a) and 6(a) representing the that is used in reconstruction, respectively. Figures 5(c)–5(e) illustrate the reconstructed surfaces using the first-order Godunov-based scheme, the high-order Godunov-based scheme, and the Lax-Friedrichs sweeping scheme, respectively. Figures 6(c)–6(e) show the corresponding reconstructed surfaces for the bottle.
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
We only evaluate the effectiveness intuitively and qualitatively. From the reconstructed results shown in Figures 5(c)–5(e), we can see that the fast eikonal solvers are more accurate than the Lax-Friedrichs sweeping scheme for mostly diffuse surface. Details of surfaces illustrated in Figures 5(c) and 5(d) are represented more vividly and clearly than in Figure 5(e). From the reconstructed results shown in Figures 6(c)–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 first-order Godunov-based scheme is the fastest and the reconstructed surface using the high-order Godunov-based scheme looks like the best result.
6. 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. The attenuation term of the light illumination will be also considered to eliminate the convex-concave ambiguity which can make the SFS problem ill-posed.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The work in this paper is supported by the Open Fund of Shaanxi Province key Laboratory of Photoelectricity Measurement and Instrument Technology, Xi’an Technological University (no. 2016SZSJ-60-1) and Xi’an Key Laboratory of Intelligent Detection and Perception (no. 201805061ZD12CG45).