Two Simple Yet Effective Strategies for Avoiding Over-Smoothing in SFS Problem

: Minimization techniques are widely used for retrieving a 3D surface starting from a single shaded image i.e., for solving the shape from shading problem. Such techniques are based on the assumption that expected surface to be retrieved coincides with the one that minimize a properly developed functional, consisting of several contributions. Among the possible contributes defining the functional, the so called “smoothness constraint” is always used since it guides the convergence of the minimization process towards a more accurate solution. Unfortunately, in areas where actually brightness changes rapidly, it also introduces an undesired over-smoothing effect. The present work proposes two simple yet effective strategies for avoiding the typical over-smoothing effect, with regards to the image regions in which this effect is particularly undesired (e.g., areas where surface details are to be preserved in the reconstruction). Tested against a set of case studies the strategies prove to outperform traditional SFS-based methods.


Introduction
One of the most used methods to retrieve the threedimensional surface of the object represented in a single image is the Shape-from-Shading (SFS) method. As widely known, SFS is an inverse problem of computer vision that, starting from the pixel by pixel analysis of the shading of an image, leads to the reconstruction of the surface of the object represented in it.
The problem, known since late '60 s (Zhang et al., 1999;Durou et al., 2008;Horn, 1970;Rindfleisch, 1966), can be presented in terms of reconstruction of the normal map of the unknown surface. A simplified formulation can be used once the following assumptions are made: • The image representing the shapes to be reconstructed is assumed to be the result of the orthogonal projection of the scene on the focal plane of the observer (i.e., perspective is absent and focal length of the observer is set at infinity) • Light beams illuminating the 2D scene are all positioned along the same, known, direction (i.e., the light source is set at infinity) • The surface is homogeneous and completely diffusing (Lambertian surface) • The represented surface does not presents hidden parts • The reference system Σ xyz , that maps the threedimensional reconstruction space, is set so that the plane П xy lies on the focal plane and the z axis is put toward the observer Under these hypotheses it is possible to formulate a relation between the surface normal N , the unknown of the reconstruction problem and the light unit-vector for each pixel of the image: Where: I(i,j) = The image, size n × m representing the shaded object whose surface is to be retrieved (input image) (i,j) = The coordinates of the generic pixel
Equation 2 regulates the irradiance of each pixel with coordinates ( , ) by modelling the reflectance of the surface as completely diffusive, i.e., without considering specular effects. Accordingly, the reconstruction problem results to be a system of irradiance equations, one for each pixel of the domain of reconstruction (generally the image without the background).
If we figure the problem graphically, Equation 2 imposes the normal N ( , ) to lay on the lateral surface of a tipped cone, whose axis coincides with the vector L and whose aperture is ( Fig. 1): Equation 3 makes more evidence the fact that the reconstruction problem is under-determined and consequently, for each pixel, the equation has infinite solutions, since there are infinite orientations for N ( , ) on the cone. This is the reason why the single image reconstruction problem using SFS, studied starting from 70 s (Zhang et al., 1999;Durou et al., 2008), has no general (closed) solution until nowadays.
In scientific literature, several approaches have been proposed to solve directly the PDEs by using propagation techniques (Kimmel and Setihian, 2004;Prados et al., 2006;Rouy and Tourin, 1992), or implementing approximation functions to both the irradiance equation and the final surface itself. However, the greater efforts in literature for solving the SFS problem are related to a range of approaches, called minimization techniques or variational methods. These methods are among the most adopted for solving the SFS problem since they prove to be extremely robust in presence of image noise or imprecise settings (e.g., guessed light direction when unknown) (Zhang et al., 1999;Durou et al., 2008;Worthington and Hancock, 1999;Huang and Smith, 2009). For this reason, such methods result applicable to a wide range of shaded images.
Variational approaches are based on the assumption that the expected surface to be retrieved, i.e., the solution that better resembles the "actual" surface represented in the shaded image, coincides with the one that minimized a properly developed functional, usually comprising the error between the (iteratively) reconstructed surface and the actual one. One of the most adopted strategies, in order to address the retrieved surface towards the expected one, is to solve the functional by imposing a number of boundary conditions. In particular, among the wide variety of boundary conditions that can be imposed to the problem, an innovative kind of them, named Morphology Based (Governi et al., 2014), proved to be effective to interactively solve the typical ambiguity between convexity and concavity on the surface. For this reason, in this work this new kind of boundary condition is used as described below.
Usually, the functional to be minimized is composed by the weighted sum of several contributions, often improperly called "constraints", each one pulling the solution towards the respect of specific requirements. The main constraints, widely used in literature, are the following: Brightness Constraint (BC), Integrability Constraint (IC) and Smoothness Constraint (SC) (Zhang et al., 1999;Durou et al., 2008;Daniel and Durou, 2000). Accordingly, in general, the functional F to be minimized is provided by the following equation: where, and are, respectively, the weights of smoothness and integrability constraints. As detailed in the following section, despite the use of a smoothness constraint in the functional is highly advisable to ensure that the minimization procedure converge to a unique solution (Zhang et al., 1999), it unfortunately also introduces possible over-smoothing effects (Worthington and Hancock, 1999). This is an undesired effect, especially for areas where, actually, brightness changes rapidly i.e., the corresponding surface is characterized by discontinuities of the shape or sharp edges (Huang and Smith, 2009;Ju et al., 2010;Chen and Dong, 2010). For this reason, the main aim of the present work is to propose two simple yet effective strategies able to avoid the typical over-smoothing effect in minimizing the SFS functional, with particular regard to the image regions in which this effect is particularly undesired (e.g., areas where brightness changes rapidly or where surface details are to be preserved in the reconstruction).

Functional Contributions
As stated in Equation 4, the functional consists of a number of constraints. The Brightness Constraint (BC), also called "variation to data" (Zhang et al., 1999;Durou et al., 2008;Daniel and Durou, 2000), is the most relevant contribution for building a functional for solving the SFS problem and, for this reason, it is always included in all the functional-based formulations. It requires that the reconstructed shape produces the same brightness as the input image at each surface point. In other words, BC is directly derived from the image irradiance and indicates the total brightness error of the retrieved surface compared with the input image (under the assumption that a constant grid of size one is considered for the image): where, R(i,j) is the estimated reflectance map and Ω is the set of all the pixel of the image (i.e., the image domain). However, the minimization of a functional comprising only this contribution would lead to an under-defined solution. This is the reason why several authors proposed a number of auxiliary constraints in order to limit the research to surfaces that satisfy not only BC but also some particular geometric or mathematical requirements. More specifically, two main kinds of constraints are, as mentioned above, available in literature: Integrability constraint and smoothness constraint.
Integrability Constraint (IC) requires the final surface to respect the principle of integrability, that limits the surface retrieval to surfaces "physically valid" (Frankot and Chellappa, 1988). From a mathematical point of view, this coincides the requirement that, for any point of the reconstruction domain, the surface height is (or needs to be) independent from the path of integration (Frankot and Chellappa, 1988;Horn, 1989;Zheng and Chellappa, 1991). Considering that, it is possible to formulate a relation between normal N and surface gradient ∇ Z as follows: Smoothness Constraint (SC) (Zhang et al., 1999;Durou et al., 2008;Governi et al., 2014;Daniel and Durou, 2000) is used to impose that the slope of the reconstructed surface changes gradually from a given pixel to its neighbourhood, so that the solution results as smooth as possible. SC is defined as follows: Consequently: The main advantage of adding this constraint into the functional is the reduction of the coarse zones on the final surface. Since minimization techniques are usually solved by imposing a number of boundary conditions (e.g., background and White points boundary conditions (Governi et al., 2013)), the SC propagates such information across the whole reconstruction domain. As already mentioned, the use of SC in the functional allows to avoid possible irregularities on the retrieved surface by limiting its roughness. Unfortunately, since it imposes a gradual change of surface normals, it might exceed in smoothing the surface, especially in areas where actually brightness rapidly changes. Consequently, the use of SC may lead to surfaces where smallest of softer details of the image are not taken into account. In other words, the over-smoothing error leads to a solution that does not satisfy the irradiance equation (Equation 1). In effect, if we observe the image obtained from the surface in the same light condition of the input image, it is possible to note that, especially in correspondence with the parts of the image in which the brightness changes rapidly, it appears blurred. Moreover, it appears more "rounded" and smoothed in correspondence with discontinuities of the shape or sharp edges, while in other zones, where the height changes more gradually, the over-smoothing effect is almost negligible (Fig. 2).

Proposed Methods
Moving from the last observation, it is possible to assert that the over-smoothing effect is concentrated in a number of zones of the surface/image, in which the slope or the brightness changes too fast.
How is it possible to overcome this problem, related to the application of the smoothness constraint?
Analysing the formulation of the smoothness constraint provided in Equation 8, it is evident that it is a sort of "chain" among the variables of the image, relating the ones relative to a given pixel to the ones relative to its neighbourhood. In the majority of the methods proposed in the scientific literature, the weight λ S of SC is set constant from pixel to pixel. This means, for instance, that a pixel on a flat surface and a pixel on the edge of a hollow equally influence the functional in terms of smoothing. Some works have taken into account this problem (Gultekin and Gokmen, 1998;Vogel et al., 2007), in particular (Vogel et al., 2007) proposed to weight the SC with a function that considers |∇I(x,y)| an edge indicator. In the same work it's proposed also the use of the evolving depth map |∇ z | as an indicator of the edges.
As previously stated, it is unthinkable to exclude the SC from the formulation of the functional, especially if dealing with the most recent (and effective) procedures, based on interactive boundary conditions (Governi et al., 2013). Consequently, the only possible way to overcome the over-smoothness is to find a method for setting the weights of SC pixel to pixel. The use of light weights for the pixels that are in proximity of a discontinuity (that may corresponds to a coarse surface) and heavy weights for the ones that lay on a smooth surface could dramatically improve the surface retrieval. Accordingly, the present work describes two main strategies, respectively called "change connectivity" and "break connectivity" for imposing variable weights to the SC in order to allow a surface reconstruction without over-smoothing. In our approach, the functional is modified as follows: Where:

Method I-Change Connectivity Strategy
This first method for imposing variable weights to the SC, named "change connectivity strategy" is based on the assumption that the "discontinuity" in the image brightness is strictly correlated to a discontinuity of the surface to be retrieved. In detail, we suppose that, in correspondence of a high variation of image brightness, there is (or, at least, there should be) a high variation in slope. Consequently, the weight of the smoothness constraint for a couple of adjacent pixels with significant difference in terms of brightness, should be less than the weight for couples of pixels with comparable brightness.
The devised approach is based on the image entropy, which is a common tool of image analysis used to evaluate the disorder of the image, which can be generated both from noise or discontinuity. As widely known the entropy H of an image is defined as follows (Shannon, 1949): where, M is the number of grey levels and p k is the probability associated with grey level k. In this study for each image pixel an entropy value is evaluated using Equation 11 applied to a v-by-v square neighbourhood (kernel): ( ) The above definition is justified by the fact that the aim here is to reduce the weight of SC in areas with greater discontinuities. Since, by definition, image entropy tends to zero if the image is uniform (flat) while it reaches its maximum value for highly disordered images, the local weight ( ) * , S i j λ evaluated for flat areas tends to 1 (i.e., high smoothness is allowed). Quite the reverse, in presence of noticeable discontinuities or edges, entropy increases thus implying the necessity of reducing the local weight to be applied to the smoothness constraint (Fig. 3). From the equation above it is evident that setting different sizes of the neighbourhood lead to different results. Moreover, the size of the neighbourhood should be chosen consistently with the size of the original image. On the basis of our experimental tests, a balanced value for v when dealing with images with maximum size 400×400 is equal to 9, since a 81 pixels kernel proves to be averagely sufficient to discriminate possible discontinuities in the image. The higher is the resolution of the original image the higher needs to be the kernel size.

Method II-Break Connectivity Strategy
The second approach proposed in this work, called "Break Connectivity Strategy" (BCS), is based on the supposition that an efficient way to "limit" the oversmoothing effect is to break the connectivity between adjacent pixels with "high" brightness gradient. By using image gradient, it is possible to isolate in the image the breaking pixels, along which propagate the fracture of the smoothness constraint. In fact, applying a gradient filter followed by a thresholding to the original shaded image it is possible to evaluate the gradient magnitude. As widely known (Gonzalez and Woods, 2008), gradient filter is a 3×3 high-pass convolution filter mainly used to detect edges in images. In particular, once the filter is applied to the generic image ( , ), a new binary image ( , ) is obtained. In such an image, edges are marked by white pixels (brightness equal to 1) while areas with low changes in slope (i.e., low brightness values) are black (brightness equal to 0). The set of white pixels in image As a consequence the resulting new formulation of smoothness constraint is the following one: In other words, pixels belonging to the domain * are characterized by a rapid change of slope or by a discontinuity with respect to a given neighborhood (e.g., abrupt uphill or downhill); for these pixels the smoothing constraint is not applied (i.e., the weight is set equal to zero).
In all the other pixels of the domain (Ω−Ω * ), the SC is equal to the one defined using a constant weight. For such a domain, even if the rate of slope is not constant from pixel to pixel, the effect of adopting a variable weight is almost unnoticeable.
Since this approach may completely separate a portion of the image from the rest of the reconstruction domain, it is crucial to pay attention in isolating only clusters of pixels in which, at least, one boundary condition is imposed. Otherwise, the solution retrieved for such partitions results incorrect, or even unfeasible and physically not valid.

Minimizing the Functional and Surface Retrieval
Once the functional is built according to Equation 9 (using one of the two strategies described above for setting the variable weights to the SC * ), the surface retrieval can be accomplished using traditional minimization techniques. In fact, the indirect minimization of the functional, aimed at evaluating a set of normals ( ) , N i j i.e., the so called "normal map", can be performed by applying literature non-linear methods. In this study, the Barzilai-Borwein non-monotonic method (Barzilai and Borwein, 1988) has been used. The unique boundary condition for constraining the minimization process is the morphology based one described in (Governi et al., 2014) and the functional is initialized using the plane normal to the light direction (Governi et al., 2014;Daniel and Durou, 2000;Governi et al., 2013). Once the normal map is evaluated, it is possible to recover the depth map of the image (i.e., the z values of the surface) using the widely known approach proposed by (Frankot and Chellappa, 1988).

Results and Discussion
Two main case studies have been carried out in order to evaluate the effectiveness of both the procedures and to figure out which strategy performs better. We used two synthetic images, i.e., two images generated directly from a CAD model representing an object characterized by a Lambertian surface properly illuminated. As a consequence-differently from the general case of SFS, in which the geometry to be retrieved is completely unknown-both height and normal map of the target surface are known a priori. This is helpful in mathematically evaluating the effectiveness of both the approaches implemented. In particular, the error between the actual (known) and the retrieved normal map is evaluated. With this aim in mind, our idea is to evaluate the angles between the true normal map and the computed one; being the normal maps composed by normal vectors, the simply scalar vector can give us the cosine of the angles: For what concerns the break connectivity strategy, given that the area of interest of this potential improvement is limited to a small bunch of pixels, it is necessary to focus on detail reproduction rather than on global relative errors like the one defined in Equation 17. For this reason, a local mean error ( ) is also used: is the domain of a 5-pixel band around the breaking pixels (composed by n b pixels).
In the following case studies, the weights relative to each constraint (SC and IC) are the following ones: λ S = 10 −1 and λ I 10 −5 . This set of weights resulted the most efficient for the reconstruction of this kind of surfaces, in this precise lighting condition.
The input image is obtained by orthogonally projecting such a 3D surface on the plane whose normal is the vector L (Fig. 4b). Since the hemisphere is obtained using a CAD software package, its actual normal map is known and therefore used as a reference for assessing the performance of the proposed methods. In particular, the two devised strategies are compared with the traditional SFS-method proposed by (Governi et al., 2013) that makes use of a constant value as a weight for smoothness constraint and the classical method proposed by Daniel Durou. It's worth to mention that Daniel-Durou's method uses a different strategy of integration from the solution given by height gradient (p, q). In fact it is used Wu and Li's method (Wu and Li, 1988) that computes the height along diagonals, then the result is used as an initial shape for Horn and Brooks' method (Horn and Brooks, 1986) which is iterative. For this reason the normals involved in the evaluation of the errors ( and ) are taken at the end of the algorithms before the integration, except for Daniel Durou algorithm where the normals are the ones given after its integration strategy. Moreover, the strategies are compared one each other.
In Fig. 5a and 5b the weight matrix Λ S and the image ( , ) for the hemisphere case study are, respectively, depicted to visually show how the weight for the SC changes in the two different strategies (CCS and BCS).
A visual comparison of retrieved surfaces is depicted in Fig. 6. In Fig. 6a a side view of the surface recovered using Durou method (Governi et al., 2014;2013) is shown. A detail of such a view is provided in Fig. 6b. In Fig. 6c and 6d are depicted, respectively, a side view of the hemisphere reconstructed using the CCS and a detail. Finally, in Fig. 6e and 6f the side view and a detail of the surface retrieved using the BCS are respectively shown. By visually comparing the three reconstructed surfaces, the best performance of the BCS leaps out, since the resulting shape is better defined and sharper around the silhouette.
The visual assessment can be confirmed by analysing the global error (Equation 17) between the actual (known) normal map and the ones obtained using, respectively, Daniel Durou algorithm, traditional method with fixed value for the weight (called "fixed lambda"), CCS and BCS methods. Such a comparison, shown in Table 1 and 2, demonstrates that the BCS method outperforms Daniel Durou (D.D. in the table), the traditional and the CCS ones. Furthermore, the performance in reconstruction using the CCS and BCS methods is, respectively, 72,27 and 76,66% better than traditional method with fixed lambda and 98,89 and 99,06% better than Daniel Durou method.       Table  3 and 4. It can be noticed that also local error in reconstruction using the proposed strategies results extremely lower with respect to traditional SFS-based reconstruction. In this case BCS proves to be 1,90% better than CCS.

Case Study 2:"Bowling Pin
The second case study consists of a figure in which a smooth surface (with a shape similar to a "bowling pin") is broken by 9 discontinuous elements (pits) i.e., discontinuities are into the reconstruction domain (and not only limited to the contour silhouette like in case study 1). In Figure 7 the CAD surface, the input image, the weight matrix and the image ( , ) for the case study are shown. It has to be noticed that for applying the BCS method only the contours of the discontinuities into the domain have been taken into account (Fig. 7d).
In Fig. 8a visual comparison between the "fixed lambda" method and the ones proposed in this work are depicted with reference to the discontinuous regions.
In Table 5 and 6 the global errors between the actual (known) normal map and the ones obtained using, respectively, Daniel Durou, fixed lambda, CCS and BCS methods are listed. In this case, it can be noticed that the global error obtained using the BCS is lower than the fixed lambda, Daniel Durou and the CCS ones, while the CCS results to be the worst one between fixed lambda and BCS methods while it's still better than Daniel Durou method. As shown in Table 7 and 8, however, the local error obtained using the CCS decreases thus demonstrating that this strategy better performs locally.  Bowling pin-global error (e g ) Fixed lambda CCS BCS 1,754481E-01 2,540573E-02 2,553353E-02 2,255123E-02 Table 6. Global error (e g ) performance increasing, see Equation 17 Bowling pin-global error (e g ) performance increasing

Case Study 3: "Golf Ball
The last case study taken into account is a sort of golf ball with 91 pits over the surface, like the ones in the bowling pin. This particular image is the biggest and has the greater number of details (e.g., holes). For these reasons the computed error is greater than the one evaluated for the previous two case studies. In Fig. 9 are shown the CAD surface, the input image, the weight matrix and the image ( , ) for this case study. It's worth to notice that for applying the BCS method in this case study, differently from the Bowling Pin case, also the background contours have been taken into account together with the contours of the discontinuities into the domain (Fig. 9d).
Also in this case study, in Table 9 and 10 it is possible to see that BCS performs better than all other methods (Fig. 10), while CCS achieves a better result than the fixed lambda and Daniel Durou methods (Fig. 11). This is evident in the local error values presented in Table 11 and 12.

Conclusion
The present paper described two simple yet effective strategies for avoiding the over-smoothing effect typically arising using the smoothness constraint for solving the SFS problem with minimization techniques. Since the smoothness constraint is used in all the literature approaches dealing with this particular topic and, to the best of our knowledge, only a few attempts in reducing the over-smoothing effect have been devised so far, the present paper could be really helpful for researchers and practitioners working in SFS field.
Test against simple case studies demonstrate the effectiveness of the two proposed strategies for surface retrieval using shaded images as input. This is particularly true when discontinuities are inside the reconstruction domain. Future work will be addressed to increase the number of test cases, with particular regard to non-synthetic and noisy images; this will allow to stress method's possible drawbacks and to conceive possible improvements. Moreover, since the BCS method seems to outperform the CCS ones, we will try to modify Equation 14 by introducing weight values in the range [0, 1] instead of constant values (0 for pixels belonging to the domain Ω * and 1 elsewhere). This could be carried out, for instance, using some outcomes of our paper (Governi et al., 2014;2013).