A New Oren–Nayar Shape-from-Shading Approach for 3D Reconstruction Using High-Order Godunov-Based Scheme

Abstract

3D shape reconstruction from images has been an important topic in the field of robot vision. Shape-From-Shading (SFS) is a classical method for determining the shape of a 3D surface from a one intensity image. The Lambertian reflectance is a fundamental assumption in conventional SFS approaches. Unfortunately, when applied to characterize the reflection attribute of the diffuse reflection, the Lambertian model is tested to be inexact. In this paper, we present a new SFS approach for 3D reconstruction of diffuse surfaces whose reflection attribute is approximated by the Oren–Nayar reflection model. The partial differential Image Irradiance Equation (IIR) is set up with this model under a single distant point light source and an orthographic camera projection whose direction coincides with the light source. Then, the IIR is converted into an eikonal equation by solving a quadratic equation that includes the 3D surface shape. The viscosity solution of the resulting eikonal equation is approximated by using the high-order Godunov-based scheme that is accelerated by means of an alternating sweeping strategy. We conduct the experiments on synthetic and real-world images, and the experimental results illustrate the effectiveness of the presented approach.

1. Introduction

3D shape reconstruction from images has been an important topic in the field of robot vision. Shape-From-shading (SFS), initiated by Horn [1,2], is a classical method for determining the shape of a 3D surface from a one intensity image. He described the SFS problem by the so-called Image Irradiance Equation (IIR), which is a nonlinear Partial Differential Equation (PDE) of the first-order, and demonstrates the relationship between the 3D surface shape and its relevant 2D image intensity variations.

There are primarily two steps when processing an SFS problem. The first step is that of setting up a partial differential IIR under certain assumptions. The assumptions are used to model the image formation process and are dependent on the reflection attribute of the surface, the light source, and the camera. The second step is that of building a numerical scheme to acquire a solution of the IIR, which is the 3D shape of the known intensity image. Since the pioneering work of Horn, the majority of the SFS works focuses on how to build a numerical scheme with the simple Lambertian model; thus, many SFS algorithms have been reported (one can see the extensive list of the references in the three surveys [3,4,5]). According to Zhang et al. [4], Durou et al. [5], and Tozza and Falcone [6], these algorithms are roughly divided into two categories: PDE-based approaches and optimization approaches. The early work of Horn [1] based on the method of characteristics and the viscosity solution theories-based methods [7,8,9,10,11] all belong to the first category. It is worth mentioning that Rouy and Tourin [7] first formulated the SFS problem as a Hamilton–Jacobi–Bellman (HJB) PDE, whose solution is obtained through the viscosity solution theories. Kimmel and Sethian [8] converted the partial differential IIR with Lambertian reflection model into an eikonal PDE. They adopted the fast-marching method [12] to compute its viscosity solution. These approaches [1,7,8] are based on an orthographic camera projection. For a perspective camera projection, Prados and Faugeras [9] associated the partial differential IIR with a Hamiltonian and approximated its viscosity solution using optimal control strategy on the basis of the work [7]. Breuß et al. [10,11] studied analytically and numerically the perspective SFS model introduced by Prados and Faugeras [9] and the related Hamilton–Jacobi (HJ) PDE. Also, they proved the convergence of a finite-difference scheme and a semi-Lagrangian scheme for the HJ PDE. The second category includes the minimization methods for the variational problem. Lee and Kuo [13] presented an iterative SFS method under the perspective camera projection. They approximated a smooth surface by the union of triangular surface patches which involved only the height variables. The shape was recovered by linearizing the reflectance map and minimizing a quadratic cost functional. Courteille et al. [14] introduced some form of regularization, by formulating a PDE-based weak perspective or perspective IIR in terms of a spline basis. We want to mention the work by Quéau et al. [15], who combined the advantages of PDE-based approaches and optimization approaches, and created a generic variational solution based on PDEs that is suitable for SFS under natural illumination and can handle a variety of scenarios for the camera and the lighting. Finally, in addition to these work, we should also mention that a latest work by Scheffler et al. [16], who gave a novel graph theoretical framework built on cycle bases. They proposed a linear program for realizing the approach to solve the disambiguation problem for SFS.

As for the modeling process of the surface reflection, few are reported in the SFS literature. Generally, the Lambertian reflectance is a fundamental assumption in conventional SFS approaches for approximating the reflection attribute of the diffuse reflection. Unfortunately, the Lambertian model has been tested to be inexact, especially for rough diffuse surfaces [17]. Samaras and Metaxas [18] adopted the Deformable Models and Oren–Nayar [17] model to construct an approach for non-Lambertian SFS and estimating the light direction. Ragheb and Hancock [19] also reported a non-Lambertian SFS with the Oren–Nayar model and presented two solutions: the analytic solution and the lookup table. Ahmed and Farag [20] gave several non-Lambertian SFS models and solved the partial differential IIR using the Lax-Friedrichs sweeping scheme [21]. Tozza and Falcone [6] proposed a unified formulation for several non-Lambertian SFS models, solved by a semi-Lagrangian approximation scheme and proved a convergence result. In these non-Lambertian SFS approaches, the distant point light source and the orthographic camera projection are adopted. Because of their simplicity, our current study is based on these models. Of course, recent research focuses on the more realistic imaging model, including a nearby point light source and perspective camera projection. Ahmed and Farag [22] used the model proposed by Prados and Faugeras [9], where the point light source is ideally attached to the camera’s projection center, and a light attenuation term 1/r² (r defines the distance between the position of the light source and the 3D surface point) has been considered in order to eliminate the ambiguity regarding convex and concave that makes the non-Lambertian SFS problem ill-posed. Moreover, Ju et al. [23] extended the work of Galliani et al. [24]: they used spherical surface parameterization to Oren–Nayar reflection model and thus could resolve an arbitrary position of the point light source. Nevertheless, the solving process needs transform the fast-marching method described in Cartesian coordinates [12] into spherical coordinates.

In this paper, we present a new Oren–Nayar SFS approach for 3D reconstruction of diffuse surfaces by extending our work [25]. The Oren–Nayar reflection model is also applied to approximate the reflection attribute of the diffuse reflection. Then, a partial differential IIR with this model is set up under a single distant point light source and an orthographic camera projection whose direction coincides with the light source. Although several works [6,20,26] have been reported in this situation, our main contribution is that we convert the IIR into a standard eikonal PDE by solving a quadratic equation which includes the 3D surface shape, after which we try to obtain the viscosity solution of the resulting eikonal equation by using the high-order Godunov-based scheme.

2. IIR for Oren–Nayar Reflectance

Under the basis that the image plane of the camera is an $x-y$ plane and the camera’s optical axis is a $z$-axis, the SFS problem is able to be expressed as reconstructing a 3D surface, $z(x,y)$, which satisfies the partial differential IIR [1]:

$$I(x,y)=R(\frac{\partial z}{\partial x},\frac{\partial z}{\partial y})=R(p(x,y),q(x,y)),$$

(1)

where $I(x,y)$ is the image irradiance that is considered to be equal to the image intensity, $p(x,y)$ and $q(x,y)$ respectively stand for the first partial derivatives of the surface $z(x,y)$ with respect to $x$ and $y$, and $R$ is the reflectance map derived from the reflection model.

Assuming that the camera performs an orthographic projection of a surface which is irradiated by a single point light source placed distantly from the surface, we can adopt the formulation as follows to describe the normal vector $\mathbf{n}$ at a 3D surface point $(x,y,z(x,y))$,

$$\mathbf{n}=[p,q,-1].$$

(2)

The reflected radiance of an ideal diffuse surface obeys Lambert’s law, which can be expressed by

$$L_r({\theta }_i)=\frac{\rho }{\pi }{I}_0\cos {\theta }_i,$$

(3)

where ${\theta }_i$, as illustrated in Figure 1, is the angle between the normal vector $\mathbf{n}$, and the incident direction of the point light source $\mathbf{L}$, $\rho$ is the albedo of diffuse surfaces that is constant in this paper, and $I_0$ is the radiant intensity of the light source. However, for a typical real diffuse surface, the Lambertian reflection model is tested to be an inexact approximation of the diffuse component of the surface reflection [17,25]. To eliminate the inexactness for diffuse reflection deriving from the assumption of the Lambertian model, Oren and Nayar [17] produced an advanced reflection model for rough diffuse surfaces where the roughness model has been used. The roughness model considers that the surface is constituted by extended symmetric V-shaped cavities, where every V-cavity consists of two planar facets which are considered to obey the Lambert’s law. The roughness of the surface is determined by a probability function of the planar facet orientations. Moreover, using a Gaussian distribution, with the reflection geometry shown in Figure 1, a formulation for the surface reflected radiance $L_r$ is addressed:

$$L_r({\varphi }_i,{\theta }_i;{\varphi }_r,{\theta }_r;\sigma )=\frac{\rho }{\pi }{I}_0\cos {\theta }_i\left(A(\sigma )+B(\sigma )\sin \alpha \tan \beta \max \left[\cos ({\varphi }_r-{\varphi }_i),0\right]\right),$$

(4)

where $\alpha =\max \left[{\theta }_i,{\theta }_r\right]$, $\beta =\min \left[{\theta }_i,{\theta }_r\right]$, $A(\sigma )=1-0.5{\sigma }^2/({\sigma }^2+0.33)$ and $B(\sigma )=0.45{\sigma }^2/({\sigma }^2+0.09)$. The parameter $\sigma$ is standed as a measure of the surface roughness, and it represents the standard deviation of the Gaussian distribution.

Considering that the direction of the point light source $\mathbf{L}$ coincides with the camera’s direction $\mathbf{V}$, we have ${\varphi }_i={\varphi }_r,{\theta }_i={\theta }_r$, and $\alpha =\beta$. Consequently, Equation (4) can be simplified to

$$L_r({\theta }_i;\sigma )=\frac{\rho }{\pi }{I}_0\left(A(\sigma )\cos {\theta }_i+B(\sigma ){\sin }^2{\theta }_i\right).$$

(5)

Generally speaking, Equation (5) can only be established under the assumptions that are a single distant point light source and an orthographic camera projection whose direction coincides with the light source or a nearby point light source attached to the projection center of the camera performing a perspective projection. If the camera’s direction does not coincide with the light source, our approach is no longer suitable, in which case one can refer to the work of [6,23]. Defining that the unit vectors of $\mathbf{L}$ and $\mathbf{V}$ are all $[0,0,-1]$ and since ${\theta }_i$ is the angle between the normal vector $\mathbf{n}$ and $\mathbf{L}$, we have

$$\begin{gathered} \cos {\theta }_i=\frac{\mathbf{L}}{\Vert \mathbf{L}\Vert }\cdot \frac{\mathbf{n}}{\Vert \mathbf{n}\Vert }=\frac{1}{\sqrt{p^2+q^2+1}}=\frac{1}{\sqrt{{\Vert \nabla z(x,y)\Vert }^2+1}}, \\ {\sin }^2{\theta }_i=1-{\cos }^2{\theta }_i=\frac{{\Vert \nabla z(x,y)\Vert }^2}{{\Vert \nabla z(x,y)\Vert }^2+1}. \end{gathered}$$

(6)

According to [1], taking Equation (6) into Equation (5), the reflectance map derived from the Oren–Nayar model is formulated as follows and now we can express the partial differential IIR for the Oren–Nayar reflectance:

$$R(p,q)=A(\sigma )\frac{1}{\sqrt{{\Vert \nabla z(x,y)\Vert }^2+1}}+B(\sigma )\frac{{\Vert \nabla z(x,y)\Vert }^2}{{\Vert \nabla z(x,y)\Vert }^2+1}=I(x,y).$$

(7)

A similar IIR for the Oren–Nayar SFS has been reported in [26] and can be solved by a semi-Lagrangian approximation scheme. Here we convert the IIR into a standard eikonal PDE which will be solved by the high-order Godunov-based scheme. Obviously, the formulation (7) is able to be characterized by a quadratic equation in the variable $\sqrt{{\Vert \nabla z(x,y)\Vert }^2+1}$. Applying the change of variable $g=\sqrt{{\Vert \nabla z(x,y)\Vert }^2+1}$, Equation (7) is formulated as

$$(I(x,y)-B(\sigma ))g^2-A(\sigma )g+B(\sigma )=0.$$

(8)

Calculating Equation (8) and satisfying $g\ge 1$, we have

$$g=\frac{2B(\sigma )}{A(\sigma )-\sqrt{A{(\sigma )}^2-4(I(x,y)-B(\sigma ))B(\sigma )}}.$$

(9)

Thus, the SFS problem (7) has been reformulated as a standard eikonal PDE:

$$\{\begin{array}{ll}\Vert \nabla z(x,y)\Vert =\sqrt{g^2-1}& \forall x\in \mathsf{\Omega },\\ z(x,y)=\varphi (x,y)& \forall x\in \partial \mathsf{\Omega },\end{array}$$

(10)

where $\varphi (x,y)$ is a boundary condition and $\mathsf{\Omega }$ is a given image domain.

3. Numerical Algorithm for Solving the Resulting Eikonal PDE

In this section, we employ the high-order Godunov-based scheme [27] that is accelerated by alternating sweeping strategy [28,29] to solve the viscosity solution of the resulting eikonal PDE (10).

3.1. High-Order Godunov-Based Scheme

Consider a uniform discretization $\{(x_m,y_n)=(ms,ns),m=1,2,\cdots ,r,n=1,2,\cdots ,c\}$ of the given image domain $\mathsf{\Omega }$ whose resolution and grid size are $r\times c$ and $s\times s$, respectively. Now our goal is to get a discrete solution $z_{m,n}=z(x_m,y_n)$ of the unknown 3D surface $z(x,y)$.

A first-order Godunov-based scheme proposed in [7,29] is employed to discretize resulting eikonal PDE (10):

$$\max {\left[\frac{z_{m,n}-z_{m-1,n}}{s},-\frac{z_{m+1,n}-z_{m,n}}{s},0\right]}^2+\max {\left[\frac{z_{m,n}-z_{m,n-1}}{s},-\frac{z_{m,n+1}-z_{m,n}}{s},0\right]}^2=F_{m,n}^2,$$

(11)

where $F_{m,n}=\sqrt{g_{m,n}^2-1}$. The substitution of $z_a=\min [z_{m-1,n},z_{m+1,n}]$, $z_b=\min [z_{m,n-1},z_{m,n+1}]$ and

$${\left\{x\right\}}^{+}=\{\begin{array}{ll}x& x\ge 0\\ 0& x<0\end{array}$$

(12)

into Equation (11), we have

$${\left[{\left\{\frac{z_{m,n}-z_a}{s}\right\}}^{+}\right]}^2+{\left[{\left\{\frac{z_{m,n}-z_b}{s}\right\}}^{+}\right]}^2=F_{m,n}^2.$$

(13)

Hence, the first-order scheme viscosity solution of Equation (13) can be obtained thusly:

$$z_{m,n}=\{\begin{array}{ll}\frac{z_a+z_b+\sqrt{2s^2{F}_{m,n}^2-{(z_a-z_b)}^2}}{2}& |z_a-z_b|<s{F}_{m,n},\\ \min [z_a,z_b]+s{F}_{m,n}& |z_a-z_b|\ge s{F}_{m,n}.\end{array}$$

(14)

To build a higher-order scheme, the $p_{m,n}$ and $q_{m,n}$ with higher-order accuracy need to be approximated. According to [27], since the third-order Weighted Essentially Non-Oscillatory (WENO) scheme has a uniform high order accuracy and is more efficient and robust than other schemes (e.g., the Essentially Non-Oscillatory (ENO) scheme), we also choose to employ the third-order WENO approximations.

$${\left\{\max \left[\frac{z_{m,n}^{new}-z_1}{s},0\right]\right\}}^2+{\left\{\max \left[\frac{z_{m,n}^{new}-z_2}{s},0\right]\right\}}^2=F_{m,n}^2,$$

(15)

where

$$\begin{array}{l}{z}_1=\min [(z_{m,n}^{old}-s{p}_{m,n}^{-}),(z_{m,n}^{old}+s{p}_{m,n}^{+})],\\ z_2=\min [(z_{m,n}^{old}-s{q}_{m,n}^{-}),(z_{m,n}^{old}+s{q}_{m,n}^{+})],\end{array}$$

(16)

with

$$\begin{array}{l}{p}_{m,n}^{-}=(1-{\alpha }_{-})\left(\frac{z_{m+1,n}-z_{m-1,n}}{2s}\right)+{\alpha }_{-}\left(\frac{3z_{m,n}-4z_{m-1,n}+z_{m-2,n}}{2s}\right),\\ p_{m,n}^{+}=(1-{\alpha }_{+})\left(\frac{z_{m+1,n}-z_{m-1,n}}{2s}\right)+{\alpha }_{+}\left(\frac{-z_{m+2,n}+4z_{m+1,n}-3z_{m,n}}{2s}\right).\end{array}$$

(17)

Here

$$\begin{array}{l}{\alpha }_{-}=\frac{1}{1+2{\beta }_{-}^2},{\beta }_{-}=\frac{\mu +{(z_{m,n}-2z_{m-1,n}+z_{m-2,n})}^2}{\mu +{(z_{m+1,n}-2z_{m,n}+z_{m-1,n})}^2},\\ {\alpha }_{+}=\frac{1}{1+2{\beta }_{+}^2},{\beta }_{+}=\frac{\mu +{(z_{m+2,n}-2z_{m+1,n}+z_{m,n})}^2}{\mu +{(z_{m+1,n}-2z_{m,n}+z_{m-1,n})}^2},\end{array}$$

(18)

where $\mu$ is a very small number keeping the denominator from getting too close to 0. In a similar way, $q_{m,n}^{-}$ and $q_{m,n}^{+}$ are specified.

Thus, the higher-order scheme viscosity solution of Equation (13) can be obtained:

$$z_{m,n}^{new}=\{\begin{array}{ll}\frac{z_1+z_2+\sqrt{2s^2{F}_{m,n}^2-{(z_1-z_2)}^2}}{2}& |z_1-z_2|<s{F}_{m,n},\\ \min [z_1,z_2]+s{F}_{m,n}& |z_1-z_2|\ge s{F}_{m,n}.\end{array}$$

(19)

3.2. Alternating Sweeping Strategy for Godunov-Based Scheme

To speed up the convergence of the high-order Godunov-based scheme, the philosophy of alternating sweeping strategy and Gauss–Seidel iterations [28,29] are used as follows. The newest available values for $z$ are adopted when we calculate the derivatives $p_{m,n}^{-}$, $p_{m,n}^{+}$, and $q_{m,n}^{-}$, $q_{m,n}^{+}$. Simultaneously, the iterations of the scheme do not only sweep in a single direction, but in the four alternating directions as follows: (a) from upper left to lower right; (b) from lower left to upper right; (c) from lower right to upper left; (d) from upper right to lower left. As can be clearly demonstrated that varied values $z_{m\pm 2,n}$, $z_{m\pm 1,n}$ and $z_{m,n\pm 2}$, $z_{m,n\pm 1}$ will be taken into account according to the current sweeping direction within the preceding cycle.

The high-order Godunov-based scheme accelerated by an alternating sweeping strategy is summarized in the following way:

(1)

Initialization: On the boundary $\partial \mathsf{\Omega }$ we set the grid points to be accurate values, i.e., $z_{m,n}^0=$ ${\varphi }_{m,n}$, which are unchanged during the iterations of the scheme. The approximated solution from the first-order Godunov-based scheme is applied as the initial guess for all other grid points.

(2)

Iterations with Alternating Sweepings Orderings: We calculate $z_{m,n}^{new}$ at $(i+1)$th iteration, updated by the rules (19) using Gauss–Seidel iterations with the following four alternating sweepings orderings:

●

$m=1:r, n=1:c$;

●

$m=r:1, n=1:c$;

●

$m=r:1, n=c:1$;

●

$m=1:r,n=c:1$.

(3)

Convergence Test: If ${\Vert z^{i+1}-z^i\Vert }_{L^1}\le \delta$, where $\delta >0$ is a known threshold value, the scheme stops and converges, or else comes back to (2). In this paper, $\delta ={10}^{-5}$ is taken.

4. Experimental Results

The effectiveness of the presented approach is assessed using several experiments on synthetic and real-world images. All the experiments were conducted on a general-purpose personal computer with an Intel Core i5-3570 CPU and 4 GB of DDR3 memory. We implemented all the approaches in Matlab, using C mex functions under Microsoft Windows operating system environments.

4.1. Experimental Results on Synthetic Images

A synthetic image of the sphere which was produced using the following formulation (20) was employed in our experiments:

$$z(x,y)=\sqrt{R^2-x^2-y^2},$$

(20)

where $(x,y)\in [-63,64]\times [-63,64]$ and $R=50$ is the sphere’s radius. Another synthetic image of the vase which was produced using the following formulation (21) was also employed:

$$z(x,y)=\sqrt{f{(y)}^2-x^2},$$

(21)

where $f(y)=-0.025(6y-1)(2y+1){(2y-1)}^2{(3y+2)}^2+0.15$ and $(x,y)\in [-0.5,0.5]\times [-0.5,0.5]$. We can map the $x$ and $y$ values to $[-63,64]$ and scale $z(x,y)$ by a factor of 128. This yields a maximum height value of approximately 36.55. Their original height maps are demonstrated in Figure 2.

The two images are generated by the model in this paper with the direction vectors of the point light source and the camera all being ${[0,0,-1]}^T$. The variant approach of Kimmel and Sethian [8] and the approach of Ahmed and Farag [20] are used to compare with our presented approach. It has been noted that the approach of Ahmed and Farag [20] applied the Lax–Friedrichs sweeping scheme to solve the Oren–Nayar SFS, and that the existing approach of Kimmel and Sethian [8] used the first-order Godunov-based scheme to solve the Lambertian SFS and therefore cannot directly solve the Oren–Nayar SFS. However, it is able to solve our partial differential IIR (10) directly. We define the first-order Godunov-based scheme that is accelerated by an alternating sweeping strategy for the Oren–Nayar SFS as the variant approach of Kimmel and Sethian [8].

The experimental results of the sphere are illustrated in Figure 3. Figure 3a shows the synthetic intensity image which is adopted for the 3D reconstruction with the parameter set: $\sigma =0.2$. Figure 3b demonstrates the reconstructed height map using the variant approach of Kimmel and Sethian [8] and the known Dirichlet boundary condition $z=\varphi$ on $\partial \mathsf{\Omega }$, while the error map with the original height map is shown in Figure 3c. Figure 3d demonstrates the reconstructed height map using the approach of Ahmed and Farag [20], while the error map with the original height map is shown in Figure 3e. Figure 3f demonstrates the reconstructed height map using our presented approach with the known Dirichlet boundary condition $z=\varphi$ on $\partial \mathsf{\Omega }$, while the error map with the original height map is shown in Figure 3g. Finally, Figure 4 illustrates the corresponding experimental results of the synthetic vase.

Figure 3 and Figure 4 show that the variant approach of Kimmel and Sethian [8], the approach of Ahmed and Farag [20], and our presented approach can all obtain satisfactory reconstructed results of the diffuse surfaces. Moreover, we can further see that our presented approach gets a higher accuracy and smaller errors than the variant approach of Kimmel and Sethian [8] and the approach of Ahmed and Farag [20].

To further demonstrate the effectiveness of our presented approach, a quantitative comparison between the variant approach of Kimmel and Sethian [8], the approach of Ahmed and Farag [20], and our presented approach is adopted by using the CPU time, the Mean Absolute Error (MAE), and the Root Mean Square Error (RMSE).

Table 1. Comparisons of approaches for the synthetic images.

Approaches	Sphere			Vase
	MAE	RMSE	CPU time	MAE	RMSE	CPU time
Variant approach of [8]	0.6903	0.8380	0.006 s	0.5464	0.6515	0.005 s
Approach of [20]	1.0637	1.4887	0.276 s	0.6852	1.0867	0.217 s
Presented approach	0.0134	0.0415	0.093 s	0.0793	0.1537	0.067 s

lists the comparisons of the three approaches for the synthetic sphere and vase images. As is directly seen, our presented approach shows better performance in the reconstructed errors than the variant approach of Kimmel and Sethian [8], as well as than the approach of Ahmed and Farag [20] in all the images which we performed. The MAE and RMSE of our proposed approach are about one-fiftieth and one-twentieth of the variant approach of Kimmel and Sethian [8] for the sphere, and about one-seventh and one-fourth for the vase, respectively. The approach of Ahmed and Farag [20] shows a slightly worse performance; maybe it is difficult to find a good estimate for the artificial viscosity term. However, our presented approach requires longer CPU time than the variant approach of Kimmel and Sethian [8], because the former adopts a complex high-order approximation scheme.

4.2. Experimental Results on Real-World Image

In order to evaluate the performance of our presented approach for real diffuse surfaces, we performed it on a real-world image and compared the reconstructed result with the existing approach of Kimmel and Sethian [8], the variant approach of Kimmel and Sethian [8], and the approach of Ahmed and Farag [20]. The image is the real-world vase used in the surveys [4,5], which is shown in Figure 5a. Figure 5b gives the ground truth referred to [5]. Figure 5c–f illustrate the reconstructed height maps using the approach of Kimmel and Sethian [8] with exact height values at each singular point, the variant approach of Kimmel and Sethian [8] with the known Dirichlet boundary condition $z=\varphi$ on $\partial \mathsf{\Omega }$, the approach of Ahmed and Farag [20], and our presented approach with the known Dirichlet boundary condition $z=\varphi$ on $\partial \mathsf{\Omega }$, respectively.

Similarly, we calculated the MAE and RMSE errors of these approaches and the comparisons results are listed in

Table 2. Comparisons of approaches for the real-world vase.

Approaches	MAE	RMSE
Approach of [8]	0.2731	0.3804
Variant approach of [8]	0.1882	0.2429
Approach of [20]	0.2648	0.3496
Presented approach	0.0771	0.0987

. From the reconstructed results illustrated in Figure 5c–f, we can see that our presented approach and the variant approach of Kimmel and Sethian [8] are more vivid than the approach of Kimmel and Sethian [8]. The same conclusion can be drawn from Table 2: our presented approach shows best performance since our approach employs the advanced reflection model for diffuse surfaces and high-order accurate numerical scheme. Recall that the approach of Kimmel and Sethian [8] assumes the Lambertian reflection model for diffuse surfaces. As previously, the approach of Ahmed and Farag [20] shows a slightly worse performance as can be seen from Figure 5e.

5. Conclusions

This paper reports a new Oren–Nayar SFS approach for 3D reconstruction of diffuse surfaces. We employed a Oren–Nayar reflection model instead of Lambertian model to approximate the reflection attribute of the diffuse reflection. The partial differential IIR for Oren–Nayar reflectance was set up with the model that the direction of the distant point light source is the same as the camera performing an orthographic projection. A standard eikonal PDE was derived from the IIR by solving a quadratic equation in the variable which includes of the 3D surface shape, and we obtained the viscosity solution of the resulting eikonal PDE using the high-order Godunov-based scheme that was accelerated by the philosophy of alternating sweeping strategy and Gauss–Seidel iterations. Experiments were conducted on synthetic and real-word images and the experimental results illustrate that our presented approach is effective and can get a higher accuracy than the first-order Godunov-based scheme and the Lax-Friedrichs scheme.

Frankly speaking, the presented approach in the paper can be only used in the situation with a single distant point light source and an orthographic camera projection whose direction coincides with the light source. Nevertheless, we believe that the idea of solving the quadratic IIR can be employed to alternately solve the SFS problem with a nearby point light source attached to the projection center of the camera performing a perspective projection. At the same time, the light attenuation term can be considered in order to eliminate the ambiguity regarding convex and concave that makes the SFS problem ill-posed. This work will be explored in the future.