Yvain Quéau
email: yvain.queau@enseeiht.fr

François Lauze
email: francois@diku.dk

Jean-Denis Durou
email: durou@irit.fr

J.-D Durou

Solving Uncalibrated Photometric Stereo using Total Variation

Keywords: 3D-reconstruction, photometric stereo, total variation, generalized bas-relief ambiguity

Estimating the shape and appearance of an object, given one or several images, is still an open and challenging research problem called 3D-reconstruction. Among the different techniques available, photometric stereo produces highly accurate results when the lighting conditions have been identified. When these conditions are unknown, the problem becomes the so-called uncalibrated photometric stereo problem, which is illposed. In this paper, we will show how total variation (TV) can be used to reduce the ambiguities of uncalibrated photometric stereo, and we will study two methods for estimating the parameters of the generalized bas-relief ambiguity. These methods will be evaluated through the 3D-reconstruction of real-world objects.

Introduction

Photometric stereo was introduced by Woodham in [START_REF] Woodham | Photometric Method for Determining Surface Orientation from Multiple Images[END_REF] in the early 80's, using the commonly adopted Lambertian model to recover both the surface shape and its albedo, given m 3 pictures of a fixed scene taken from the same viewpoint but under different calibrated illumination conditions. However, calibrating the light sources is not a trivial task, and it has been shown that slight errors in this calibration could result in major reconstruction errors [START_REF] Horovitz | Depth from gradient fields and control points: Bias correction in photometric stereo[END_REF]. Thus, uncalibrated photometric stereo seems to be a promising alternative, as the light vectors are estimated during the reconstruction. Unfortunately, this problem is ill-posed. If integrability is imposed, the ambiguity is reduced to a set of 3-parameter linear transformations called generalized bas-relief transformations [START_REF] Belhumeur | The Bas-Relief Ambiguity[END_REF].

The estimation of these parameters is the subject of this paper: we extend the preliminary results presented in [START_REF] Quéau | Solving the Uncalibrated Photometric Stereo Problem Using Total Variation[END_REF] and show how TV-regularization of the M field, which encodes the albedo and the normals, or of the depth function u, considerably reduces the ambiguities: the generalized bas-relief ambiguity is reduced to the classical bas-relief ambiguity, which is easier to solve. We demonstrate that the proposed methods can be written as the minimization of simple convex functions, and are thus much faster than most previous work involving either pixel identification [START_REF] Favaro | A Closed-Form Solution to Uncalibrated Photometric Stereo via Diffuse Maxima[END_REF][START_REF] Shi | Self-calibrating Photometric Stereo[END_REF] or nonconvex minimization [START_REF] Alldrin | Resolving the Generalized Bas-relief Ambiguity by Entropy Minimization[END_REF]. Despite the absence of statistical justification on the distributions of the albedo, the depth and the normal field, the results we obtain a posteriori justify the choice of using total variation.

The rest of the paper is organized as follows: subsequently to the recall of the equations of Lambertian photometric stereo in Section 2, we empirically justify in Section 3 the choice of using TV for solving the uncalibrated photometric stereo problem. In Section 4, we prove that minimizing the total variation of the M field or that of the depth u reduces the generalized bas-relief ambiguity to a simple bas-relief ambiguity. The methods are eventually evaluated in Section 51 .

Background

Calibrated Photometric Stereo

Photometric Stereo (PS) was introduced in [START_REF] Woodham | Photometric Method for Determining Surface Orientation from Multiple Images[END_REF] to deal with 3D-reconstruction, allowing the user to recover both the normal field associated to a surface and the albedo. In this technique, m images of a scene are taken from the same viewpoint but under variable lighting conditions. An example of such images, with m = 3, is shown in Fig. 1. PS is an extension of the shape-from-shading problem [START_REF] Horn | Shape from Shading[END_REF], which is known for being ill-posed [START_REF] Oliensis | Uniqueness in shape from shading[END_REF]. The use of additional images with different lighting conditions allows to solving the ambiguities of shape-fromshading. Instead of a unique known light source, supposed to be located at infinite or finite distance (a recent shape-from-shading considered arbitrary location [START_REF] Galliani | Generalised Perspective Shape from Shading in Spherical Coordinates[END_REF]), dealing with several known sources allows obtaining unambiguous, accurate and fast 3D-reconstructions. As the normal vectors and the albedo values are estimated at every pixel, dense depth and albedo maps are recovered, unlike using stereoscopy. To obtain a full 3D-model rather than a depth map, PS can be coupled with multi-view techniques, through the fusion of depth maps [START_REF] Chang | Robust Estimation and Organization of Shapes from Multiple Cues 14[END_REF] or the combination with, e.g., structure-frommotion [START_REF] Zhou | Multi-view Photometric Stereo with Spatially Varying Isotropic Materials[END_REF] or shape-from-silhouettes [START_REF] Hernández | Multiview Photometric Stereo[END_REF].

Assuming orthographic projection (the problem was also studied in the perspective case [START_REF] Mecca | Perspective Photometric Stereo with Shadows[END_REF]), the surface to be reconstructed is represented by a Monge patch (x, y, u(x, y)) ⊤ over a domain Ω ⊂ R 2 , where the depth function u(x, y) is assumed to be twice differentiable, at least piecewise. The vector field of unit outward normals is given by N (x, y) = 1

∇u 2 + 1   -∂ x u -∂ y u 1   (1)
Let ρ(x, y) ∈ [0, 1] be the albedo map, and M (x, y) = ρ(x, y)N (x, y), so that ρ(x, y) = ||M (x, y)||. Given a 2 http://www.ece.ncsu.edu/imaging/Archives/ImageDatabase/ directional light source S i = (S i x , S i y , S i z ) ⊤ ∈ R 3 , the graylevel I i (x, y) obtained at point (x, y) of the image plane for matte objects is given by Lambert's law I i (x, y) = ρ(x, y)N (x, y) ⊤ S i = M (x, y) ⊤ S i

(2)

(shadows are ignored).

In order to recover M (and thus ρ and N ) without ambiguity, at least 3 images associated with 3 noncoplanar light vectors are mandatory. In the sequel the m 3 images will be denoted I 1 , . . . , I m , and the corresponding light vectors S 1 , . . . , S m . It will also be assumed that ∀(x, y), ρ(x, y) > 0. Defining the observation vector I(x, y) = (I 1 (x, y), . . . , I m (x, y)), the light matrix S = (S 1 , . . . , S m ), and omitting the dependencies in (x, y), the following equality holds:

I = M ⊤ S (3)
This is a system of linear equations which can be solved in a least-squares sense using the Moore-Penrose pseudoinverse S + of S:

M ⊤ = IS + , ρ = M , N = M ρ (4)
(all through this paper, a hat above a letter will indicate that this is an estimation).

The normal vectors being unitary, they have only two degrees of freedom, so the problem is often equivalently formulated as the estimation of ρ and ∇u = (p, q) ⊤ , with

     p = - N x N z q = - N y N z (5)
The first line of Fig. 2 shows the albedo and the depth gradient estimated from the images of Fig. 1.

Although direct depth estimation may be carried out [START_REF] Mecca | Perspective Photometric Stereo with Shadows[END_REF], an additional step is usually needed to "integrate" the normal field into a depth map (see Fig. 2-d). This step, which can be problematic, will not be presented in this paper. It usually results in estimating the depth u as the least-squares solution of the equation ∇u = ( p, q) ⊤ [START_REF] Durou | Integration of a normal field without boundary condition[END_REF][START_REF] Frankot | A Method for Enforcing Integrability in Shape from Shading Algorithms[END_REF][START_REF] Harker | Least squares surface reconstruction from measured gradient fields[END_REF][START_REF] Horn | The variational approach to shape from shading[END_REF][START_REF] Simchony | Direct analytical methods for solving Poisson equations in computer vision problems[END_REF]. All the depth maps presented in this paper are obtained using the same solver [START_REF] Simchony | Direct analytical methods for solving Poisson equations in computer vision problems[END_REF].

About the Assumptions

The usual assumption about the number of images (i.e. the number of light sources) is that m 3 and the light vectors should be non-coplanar. Some attempts have been made to solve the problem with m = 2 [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF][START_REF] Mecca | Unambiguous photometric stereo using two images[END_REF][START_REF] Onn | Integrability Disambiguates Surface Recovery in Two-Image Photometric Stereo[END_REF], and an interesting recent application of the coplanar configuration is 3D-reconstruction from an outdoor webcam (since the Sun moves within a plane) [START_REF] Abrams | Heliometric Stereo: Shape from Sun Position[END_REF][START_REF] Ackermann | Photometric Stereo for Outdoor Webcams[END_REF]. In the present paper the traditional case of m 3 noncoplanar light vectors is studied.

In addition, although the m > 3 case improves the robustness of PS thanks to the least-squares estimation (4), a lot of physical phenomena like specularities, self-shadows, cast-shadows or inter-reflections are not represented in Lambert's law (2): such phenomena can be treated as outliers in a robust estimation process like Expectation-Maximization [START_REF] Verbiest | Photometric stereo with coherent outlier handling and confidence estimation[END_REF]. Specular reflections have also been successfully handled by adapting the photometric stereo technique to non-purely diffuse models such as that of Torrance and Sparrow [START_REF] Georghiades | Incorporating the Torrance and Sparrow model of reflectance in uncalibrated photometric stereo[END_REF], considering that the luminance is obtained by summing a diffuse and a specular component.

The problem of self-shadows, specularities, and more generally arbitrary reflectances can also be avoided by using a reference object [START_REF] Hertzmann | Example-Based Photometric Stereo: Shape Reconstruction with General, Varying BRDFs[END_REF] whose geometry and reflectance are known. The case of cast-shadows is more complicated, since they are impossible to locally model, however they can be located and rejected from the estimation process [START_REF] Barsky | The 4-source photometric stereo technique for three-dimensional surfaces in the presence of highlights and shadows[END_REF]. This has been done efficiently using the popular Graph Cuts technique [START_REF] Chandraker | Reflections on the Generalized Bas-Relief Ambiguity[END_REF].

Another approach to deal with outliers is to preprocess the data: under the Lambertian assumption, any set of n pixels of the observed images forms a n × m matrix which should be of rank 3. However, this is not the case when dealing with real images, because of noise, quantization and outliers. Thus, Wu et al. propose in [START_REF] Wu | Robust Photometric Stereo via Low-Rank Matrix Completion and Recovery[END_REF] an algorithm to correct the input, via low-rank approximation. It should also be noted that specularities can be removed from RGB images by an appropriate change of colour space [START_REF] Mallick | Beyond lambert: Reconstructing specular surfaces using color[END_REF]. In this work, the outliers are ignored, i.e. it is assumed that we are dealing either with synthetic data or with preprocessed real data.

Uncalibrated Photometric Stereo

In Section 4, we will focus on the problem with m 3 unknown non-coplanar distant light sources, which is the so-called uncalibrated photometric stereo problem. It this case, it has been shown [START_REF] Epstein | Learning and Recognizing Objects Using Illumination Subspaces[END_REF][START_REF] Hayakawa | Photometric stereo under a light-source with arbitrary motion[END_REF] that the vector field M defined as the surface normal multiplied by the albedo can only be determined up to a 3 × 3 linear transformation A ∈ GL(3), since

M ⊤ S = A ⊤ M ⊤ A -1 S (6)
However, since the depth function u(x, y) is assumed to be twice differentiable, Schwarz' integrability constraint ∂ 2 u/∂ x ∂ y = ∂ 2 u/∂ y ∂ x should be satisfied. This condition is equivalent to equating to zero the third component of curl (N/N z ) = curl (M/M z ), denoted by curl M :

curl M = ∂ ∂y M x M z - ∂ ∂x M y M z = 0 (7)
Imposing curl M = 0 as a hard constraint is not possible, as Eq. ( 7) does not hold for points located at discontinuities, which are present in real-world objects. Thus, Yuille and Snow proposed in [49] to solve Eq. ( 7) in a least-squares sense, while respecting Lambert's law (3) (still in a least-squares sense). We reformulate this method as the estimation of a field M and a light matrix S holding both constraints:

M ⊤ S -I 2 = min curl M 2 = min ( 8
)
where the notation = min signifies that the parameters minimize the left hand of the equation. Adapting Hayakawa's method [START_REF] Hayakawa | Photometric stereo under a light-source with arbitrary motion[END_REF], Yuille and Snow give in [49] a closed-form solution holding these constraints, which will be referred to as (M 0 , S 0 ) in the following. Not all transformations A ∈ GL(3) preserve the integrability of the field A ⊤ M : it has been shown in [START_REF] Belhumeur | The Bas-Relief Ambiguity[END_REF]49] that imposing the integrability constraint reduces the set of transformations A to the set of 3-parameters transformations A = G with G of the form

G(µ, ν, λ) =   1 0 0 0 1 0 µ ν λ   , G(µ, ν, λ) -1 = 1 λ   λ 0 0 0 λ 0 -µ -ν 1   (9)
where (µ, ν, λ) ∈ R 2 × R\{0}. This corresponds to the generalized bas-relief (GBR) ambiguity [START_REF] Belhumeur | The Bas-Relief Ambiguity[END_REF]. As the sign of λ is directly linked to the convexity of the object, and as the concave/convex ambiguity cannot be solved without additional information on the object, it is usually assumed that λ > 0. Thus, problem ( 8) is still an ill-posed problem, as one can apply any GBR transformation without either changing the value of the data term (this comes from ( 6)), or breaking the integrability condition, since it is shown in [START_REF] Belhumeur | The Bas-Relief Ambiguity[END_REF]49]

that curl M 2 = min ⇐⇒ curl(G ⊤ M ) 2 = min (10)
In terms of depth, the GBR transforms u(x, y) into ū(x, y) = u(x, y) -µx -νy λ (11) Fig. 3 shows the effects of a GBR.

Research on the resolution of this GBR ambiguity is prolific. It can be assumed, quite reasonably, that real-world objects do not strictly follow the Lambertian model, so as to exploit the presence of outliers like specularities [14] or inter-reflections [START_REF] Chandraker | Reflections on the Generalized Bas-Relief Ambiguity[END_REF]. Local properties of the shape or the reflectance can also be used: for example, critical points in an image (the brightest pixels) can reveal the direction of the light vector [START_REF] Favaro | A Closed-Form Solution to Uncalibrated Photometric Stereo via Diffuse Maxima[END_REF], and intensity profiles allow identifying pixels which have the same normal but different albedos [START_REF] Shi | Self-calibrating Photometric Stereo[END_REF]. We argue that considering global properties is more general, and might be easier to exploit as it does not rely on the identification of critical points.

Quite surprisingly, very few attemps to exploit global properties have been made, apart from the use of prior on shapes [START_REF] Georghiades | From Few to Many: Illumination Cone Models for Face Recognition under Variable Lighting and Pose[END_REF] or on the albedo distribution: Alldrin et al. propose in [START_REF] Alldrin | Resolving the Generalized Bas-relief Ambiguity by Entropy Minimization[END_REF] to estimate the GBR parameters which minimize the entropy of the albedo distribution, since the albedo is spread by the GBR. Their motivation is to favour materials which are "homogeneous", i.e. made up of a small value of components. Consistently with our formulation (8), we rewrite their method The scales are the same as in Fig. 2, so as to illustrate the drifts in all the three directions. It is to be noted that the albedo (a) is spread, compared to that of Fig. 2-a: this led to the Minimum-Entropy algorithm [START_REF] Alldrin | Resolving the Generalized Bas-relief Ambiguity by Entropy Minimization[END_REF].

as the following pair of problems, which must be solved sequentially:

       M 0 ⊤ S 0 -I 2 = min curl M 0 2 = min ( µ, ν, λ) = argmin µ,ν,λ E( G(µ, ν, λ) ⊤ M 0 ) ( 12
)
where E( G ⊤ (µ, ν, λ)M 0 ) stands for the entropy of the transformed albedo distribution.

TV-Regularized Models in Computer Vision

The total variation of a function is a widely used measure for regularity. For an almost everywhere differentiable function f : Ω ⊂ R p → R, it can be written

TV(f ) = Ω ∇f (x) dx (13)
and extends to the class of so-called functions of bounded variations [START_REF] Aubert | Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations[END_REF].

In the sequel we take p = 2 as we deal with images. Total variation was introduced to the imaging community by the work of Rudin, Osher and Fatemi (ROF) [START_REF] Rudin | Nonlinear Total Variation Based Noise Removal Algorithms[END_REF] for image denosing: given a (noisy) graylevel image I 0 = I + η, with η a Gaussian white noise of known variance σ 2 , the ROF model finds an image I which solves I -I 0 2 ≤ σ 2 on Ω, and with minimal total variation T V ( I). The problem was shown in [START_REF] Chambolle | Image recovery via total variation minimization and related problems[END_REF] to be equivalent to the minimization of the functional

E(I) = I -I 0 2 + γ σ T V (I)
for a certain weight γ σ . Efficient algorithms have been developed in the last decade for solving such problems, using for instance primal-dual schemes [START_REF] Chambolle | An Algorithm for Total Variation Minimization and Applications[END_REF]. In this paper, since both the total variation of the depth u and that of the R 2 → R 3 vector field M will be considered, vectorial TV should also be defined. When f takes its values in R q with q > 1, several definitions can be considered (see [START_REF] Bresson | Fast Dual Minimization of the Vectorial Total Variation Norm and Applications to Color Image Processing[END_REF] for some discussion). The most frequently used definition is probably

TV(f ) = Ω J(f (x)) F dx (14)
where J(f (x)) F is the Frobenius norm of the Jacobian matrix of f = (f 1 . . . f q ) ⊤ at point x. Another definition is:

TV(f ) = q i=1 Ω ∇f i (x) dx = q i=1 TV(f i ) ( 15
)
The difference between both these definitions may be important in the denoising context, as the first introduces a coupling between the different channels, and is invariant by rotation of these channels. However, in our context, we are not concerned by this difference, as the ordering of the three components M x , M y and M z is important for the integrability constraint [START_REF] Belhumeur | The Bas-Relief Ambiguity[END_REF], and as each parameter of a GBR affects a single component of M . Therefore, we will consider the definition [START_REF] Durou | Integration of a normal field without boundary condition[END_REF], which is easier to manipulate and produces better results. In this paper, unlike the preliminary version [START_REF] Quéau | Solving the Uncalibrated Photometric Stereo Problem Using Total Variation[END_REF], the total variation of the M field will thus refer to

TV(M ) = TV(M x ) + TV(M y ) + TV(M z ) (16)
In Section 4, we will show that minimizing the total variation of either the M field or the depth function u reduces the ambiguities of uncalibrated photometric stereo. As we will see, in this context, the minimization of TV is quite simple as it depends only on the GBR parameters. The resulting problems being convex, standard convex optimization methods such as quasi-Newton ones can be used: this offers considerably faster methods, compared to state-of-the-art approaches.

An Empirical Evidence for using TV

Let us now introduce our motivations for using TV in the context of uncalibrated photometric stereo. As it will be illustrated on some examples, a GBR affects the variations of both the albedo ρ, the depth u and its gradient ∇u, which is linked to the normals by [START_REF] Abrams | Heliometric Stereo: Shape from Sun Position[END_REF]. Moreover, the optimal values of µ and ν seem to correspond to a minimum of the total variation.

Firstly, let us state that this work is inspired by the method of Alldrin et al., who show in [START_REF] Alldrin | Resolving the Generalized Bas-relief Ambiguity by Entropy Minimization[END_REF] that the GBR parameters could be estimated by minimizing the entropy of the albedo distribution. However, the entropy being anything but convex, the minimization has to be global: Alldrin et al. use a bruteforce discrete search on a user-defined interval, which results in a very long optimization process. Furthermore, this entropy does not consider the spatial variation of the albedo, as shown in Fig. 4. When looking for "homogeneous" zones, one would expect that similar albedo pixels would be close to each other. This property can be obtained by using the total variation of the albedo, as one would expect that locally, the variations of albedo should be slight. Also, total variation having better differential properties than entropy, we expect considerably shorter computation times since convex optimization tools can be used instead of bruteforce algorithms.

Since the GBR applies not only to the albedo but to the whole M field, not only the variations of the albedo ρ are affected but also those of the depth u and of its first order derivatives ∇u = (p, q) ⊤ . We show in Fig. 5 the effect of a GBR on these variations. We plotted the values of the ground truth and the GBR-transformed values of (e) ρ, (f) u, (g) p and (h) q along the y = 100 straight line, with µ = 1/2, ν = 1 and λ = 1/3. Obviously, the GBR-transformed estimations have higher total variations. Thus, to solve this GBR ambiguity, it seems reasonable to consider not only ρ, but also u and ∇u. As ρ, ∇u and M are linked by

M = ρ ∇u 2 + 1 -∇u 1 ( 17
)
we can also directly consider the M field.

In Fig. 6, we report the impact of µ and ν on the total variations of ρ, u and M estimated from the images of Beethoven's bust. To do so, we consider as the ground truth the solution of the calibrated photometric stereo problem and apply GBR transformations to this ground truth, with λ fixed to 1: the ground truth values thus correspond to (µ, ν) = (0, 0). To independently study the impact of each parameter, one test per parameter is executed while fixing the other parameters to the ground truth value. This test tends to indicate that minimizing any of these total variations may solve the problem in µ and ν. The case of λ is different: the effect of λ on these TV, with µ and ν fixed to the ground truth values, is illustrated in Fig. 7. TV(u) is minimized for λ = +∞ and TV(M ) for λ = 0 (these observations are proven in Section 4), however these values correspond to degenerated depth maps (see Eq. ( 11)). Thus, only TV(ρ) seems to have a non-degenerated minimum. However, in general, estimating λ by minimizing TV(ρ) gives quite disappointing results: in fact the previous test is biased by the obviously uniform albedo of Beethoven's bust (which is made of plaster). Using the dataset Doll from [START_REF] Alldrin | Resolving the Generalized Bas-relief Ambiguity by Entropy Minimization[END_REF], we obtain the plot of Fig. 8 for TV(ρ) as a function of λ, with µ and ν fixed to the ground truth values. Thus, to estimate λ, other methods should be considered. However, if µ and ν are already estimated, we will see in Section 4.3 that this estimation is equivalent to the resolution of the classical bas-relief ambiguity, which is considerably easier to solve than the GBR. As a consequence, very simple solutions do exist. For instance, we will show that the constant light magnitude constraint [START_REF] Hayakawa | Photometric stereo under a light-source with arbitrary motion[END_REF]49] can be written as a simple linear system. Of course, any method of resolution of the GBR [START_REF] Alldrin | Resolving the Generalized Bas-relief Ambiguity by Entropy Minimization[END_REF][START_REF] Favaro | A Closed-Form Solution to Uncalibrated Photometric Stereo via Diffuse Maxima[END_REF][START_REF] Shi | Self-calibrating Photometric Stereo[END_REF] would also be simplified by the knowledge of µ and ν (see for instance the TV-u / ME method in Section 5). Because of the link (17) between ρ, u and M , we restrict ourselves in this paper to the study of TV(M ) and TV(u). Moreover, the following proposition shows that minimizing TV(M ) is linked (although not equivalent) to simultaneously minimizing TV(ρ) and TV(∇u):

Proposition 1 TV(M ) 3TV(ρ) + 4TV(∇u)
Proof Denoting ∇u = (p, q) ⊤ and r = ∇u , we have

M x = -ρp √ r 2 + 1 , M y = -ρq √ r 2 + 1 , M z = ρ √ r 2 + 1 ( 18
)
The gradient of M x is given by

∇M x = - ρ∇p √ r 2 + 1 - p∇ρ √ r 2 + 1 + ρpr∇r (r 2 + 1) 3/2 (19) (20) Since 1 √ r 2 +1 1, p √ r 2 +1
1 and r r 2 +1

1, we get

∇M x ρ ∇p + ∇ρ + ρ ∇r (21)
We obtain in the same way

∇M y ρ ∇q + ∇ρ + ρ ∇r (22)
The gradient of M z is given by

∇M z = ∇ρ √ r 2 + 1 - ρr∇r (r 2 + 1) 3/2 (23)
Thus we have

∇M z ∇ρ + ρ ∇r (24)
Summing ( 21), ( 22) and ( 24) gives c=x,y,z

∇M c 3 ∇ρ + 3ρ ∇r + ρ( ∇p + ∇q ) (25)
Finally, since ∇r = p∇p+q∇q r and max(|p|, |q|) r,

∇r ∇p + ∇q (26)
and thus c=x,y,z

∇M c 3 ∇ρ + 4ρ( ∇p + ∇q ) (27)
Integrating this inequality over the reconstruction domain Ω leads to:

TV(M ) 3 Ω ∇ρ dx + 4 Ω ρ( ∇p + ∇q )dx (28)
Since ρ ∈ [0, 1], using the definition (13) of TV(ρ) and the definition [START_REF] Durou | Integration of a normal field without boundary condition[END_REF] of TV(∇u) leads to TV(M ) 3 TV(ρ) + 4 TV(∇u) [START_REF] Horn | Shape from Shading[END_REF] which is the announced result. ⊓ ⊔ Remark 1 In Eq. ( 28), Ω ρ( ∇p + ∇q )dx is the ρweighted total variation of ∇u and when minimizing it, it allows for "relaxing" the minimization of TV(∇u) where the albedo is low, i.e. where the material is dark. Since it is obvious that an albedo equal to zero induces ill-posedness, this "relaxation" allows us not to consider areas which would induce errors in the reconstruction.

Finally, let us state that, though the proposed limited empirical evidence should be completed by a statistical study of the involved distributions, the results shown in Section 5 will a posteriori justify the TV approach, as the accuracy of the estimations of µ and ν reach state-of-the-art results, and the computation times of the proposed methods are considerably shorter than those of other approaches, thanks to the convexity of the objective functions that are discussed in the next section.

Solving the GBR using TV

In this section, we show how TV-regularization can solve the uncalibrated photometric stereo problem. Our aim is to estimate a solution ( M , S) of ( 8) having minimal total variation. Given a solution (M 0 , S 0 ) of ( 8), this problem can be defined as the estimation of the GBR parameters µ, ν and λ: the M field is eventually given by M = G( µ, ν, λ) ⊤ M 0 .

The Proposed Model

In order to minimize the TV of the M field, different approaches can be considered. The first, inspired by denoising models, would be to replace (8) by:

   ( M , S) = argmin M,S M ⊤ S -I 2 + γ TV(M ) s.t. curl M 2 = min ( 30
)
where γ > 0 is a user-defined hyper-parameter. However, this is not the best solution: such a model would be useful in order to smooth some artifacts due to the presence of outliers, however, as already stated, in this paper we are not concerned about such outliers as we assume the data is preprocessed. Therefore, a weight between the data term M ⊤ S -I 2 and the regularization term TV(M ) is not required, since we would like to keep the data term minimal 3 .

In the context of the resolution of the GBR ambiguity, we propose to solve the following problem, which seems more relevant:

       M = argmin M TV(M ) s.t. M ⊤ S -I 2 = min curl M 2 = min (31)
To satisfy the first constraint in [START_REF] Koenderink | The generic bilinear calibration-estimation problem[END_REF], orthogonal projection of a local least-squares estimation on the space of curl-free vector field might be considered, as such a projection can be achieved efficiently using, for instance, Fourier transform [START_REF] Frankot | A Method for Enforcing Integrability in Shape from Shading Algorithms[END_REF][START_REF] Simchony | Direct analytical methods for solving Poisson equations in computer vision problems[END_REF]. Unfortunately, since the "curl" considered here (defined by Eq. ( 7)) is nonlinear, such an approach is not as simple to put in place as it may seem. Another possibility would be to introduce two Lagrange multipliers associated with the constraints. As the second constraint is nonlinear and as there are several unknowns, proximal methods [START_REF] Attouch | Alternating Proximal Algorithms for Weakly Coupled Convex Minimization Problems. Applications to Dynamical Games and PDE's[END_REF][START_REF] Chambolle | A first-order primaldual algorithm for convex problems with applications to imaging[END_REF] could probably be used. Yet, this would require an iterative scheme, which would result in high computation times.

Moreover, we did not use the fact that we do know a closed form solution of (8), and the analytical form of all acceptable transformations (the GBR matrices). We propose a simple solution to [START_REF] Koenderink | The generic bilinear calibration-estimation problem[END_REF], based on the solution of Yuille and Snow [49] and the estimation of the GBR parameters. We use the following proposition:

Proposition 2 Problem (31) is strictly equivalent to ( M , S) = G( µ, ν, λ) ⊤ M 0 , G( µ, ν, λ) -1 S 0 (32)
where (M 0 , S 0 ) is the solution of Yuille and Snow, and

   ( µ, ν, λ) = argmin µ,ν,λ TV(G(µ, ν, λ) ⊤ M 0 ) s.t. λ > 0 ( 33
)
Proof Let ( M , S) be given by ( 32) and [START_REF] Lee | Acquiring Linear Subspaces for Face Recognition under Variable Lighting[END_REF]. A method for estimating µ, ν and λ is proposed in the following section. All the GBR matrices being invertible, we have

M 0 ⊤ S 0 -I 2 = (G( µ, ν, λ) ⊤ M 0 ) ⊤ (G( µ, ν, λ) -1 S 0 )-I 2 .
As, by construction, M 0 ⊤ S 0 -I 2 = min, it follows that M ⊤ S -I 2 = min. By construction also, M 0 is such that curl M 0 2 = min. From (10), we deduce curl M 2 = min. Both constraints are thus respected, thanks to the properties of the GBR. Finally, the GBR parameters being chosen so as to minimize TV(G(µ, ν, λ) ⊤ M 0 ), M is a solution of [START_REF] Koenderink | The generic bilinear calibration-estimation problem[END_REF].

Now, let M be a solution of [START_REF] Koenderink | The generic bilinear calibration-estimation problem[END_REF], and S the associated light matrix. ( M , S) satisfies the two constraints, as does (M 0 , S 0 ). As the only transformations which preserve these constraints are the GBR [START_REF] Belhumeur | The Bas-Relief Ambiguity[END_REF]49], there necessarily exists a GBR matrix G( µ, ν, λ) such that M = G( µ, ν, λ) ⊤ M 0 and S = G( µ, ν, λ) -1 S 0 , i.e. ( M , S) can be expressed under the form [START_REF] Kozera | Existence and uniqueness in photometric stereo[END_REF].

⊓ ⊔

Let us now explain how to solve (33).

Estimating µ and ν

Firstly, for computational issues, we propose to replace the constraint λ > 0 by λ ǫ, with ǫ > 0. We will discuss the influence of ǫ in Section 4.3, and propose solutions for choosing it (in the preliminary results [START_REF] Quéau | Solving the Uncalibrated Photometric Stereo Problem Using Total Variation[END_REF], we manually chose ǫ). Now, remark that, from our definition (15) of TV and the form of a GBR, the objective function reads:

E(µ, ν, λ) = TV(G(µ, ν, λ) ⊤ M 0 ) = TV(g 1 (µ)M 0 ) + TV(g 2 (ν)M 0 ) + TV(g 3 (λ)M 0 ) ( 34
)
where g 1 (µ) = (1, 0, µ), g 2 (ν) = (0, 1, ν) and g 3 (λ) = (0, 0, λ) denote the three lines of G(µ, ν, λ) ⊤ .

Writing TV(g 3 (λ)M 0 ) = Ω ∇(λM 0 z ) dx, one obtains ∂ λ E = λ TV(M 0 z ), and thus TV(M ) is linear in λ (cf. Fig. 7-c), and the solution in λ is given by:

λ = ǫ (35)
This explains why we replaced the constraint λ > 0 by λ ǫ, as λ = 0 would provide a degenerated depth function ū(x, y) (see [START_REF] Chambolle | A first-order primaldual algorithm for convex problems with applications to imaging[END_REF]). Thus, the choice of ǫ is important for the estimation of λ, but we can state the following:

Proposition 3 The estimations of µ and ν are totally independent from that of λ, thus from ǫ.

Proof Denoting by E(µ, ν, λ) = TV(G(µ, ν, λ) ⊤ M 0 ) = TV(M 0
x + µM 0 z ) + TV(M 0 y + νM 0 z ) + TV(λM 0 z ), its derivatives in µ and ν read:

       ∂ µ E = Ω µ ∇M 0 z 2 + (∇M 0 x ) ⊤ (∇M 0 z ) ∇(M 0 x + µM 0 z ) dx ∂ ν E = Ω ν ∇M 0 z 2 + (∇M 0 y ) ⊤ (∇M 0 z ) ∇(M 0 y + νM 0 z ) dx (36)
Since λ does not appear in Eq. ( 36), this proves the proposition. It should be stated that this comes from the choice of the definition (15) of TV: λ would remain in the equations if the definition (14) were used.

⊓ ⊔

Remark 2 Differentiating Eq. ( 36) w.r.t. µ and ν proves the convexity of TV(M ) w.r.t. µ and ν.

As will be shown in Section 4.3, the choice of ǫ actually corresponds to the resolution of the classical bas-relief ambiguity, which can be solved quite simply. Hence, the proposed TV formulation reduces the GBR to the bas-relief ambiguity. As already stated by Belhumeur et al. in [START_REF] Belhumeur | The Bas-Relief Ambiguity[END_REF], this is consistent with the work of Fan and Wolff [START_REF] Fan | Surface curvature and shape reconstruction from unknown multiple illumination and integrability[END_REF], who show that the Hessian of u can only be recovered up to a scaling factor (which here corresponds to λ). In their work, two other parameters are required (the boundary conditions of a differential equation): these parameters correspond to the remaining parameters µ and ν.

At this point, let us focus on the estimations of µ and ν. They involve the minimization of a convex function which depends only on two parameters µ and ν. Given the expression (36) of the gradient ∇E, we solve for µ and ν in problem (33) using a quasi-Newton method. In our implementation, we use the fminunc Matlab function. We can accelerate the optimization by choosing an adequate initialization, such as the solution of the Tikhonov regularization, given by

     µ 0 = argmin µ Ω ∇(M 0 x + µM 0 z ) 2 dx ν 0 = argmin ν Ω ∇(M 0 y + νM 0 z ) 2 dx ( 37
)
which admits the closed-form solution

         0 = -Ω (∇M 0 x ) ⊤ (∇M 0 z )dx Ω ∇M 0 z 2 dx ν 0 = -Ω (∇M 0 y ) ⊤ (∇M 0 z )dx Ω ∇M 0 z 2 dx (38)

Estimating λ

According to [START_REF] Mecca | Unambiguous photometric stereo using two images[END_REF], the solution in λ of problem ( 33) is λ = ǫ. Similarly to a hyper-parameter, ǫ can be manually tuned by the user so as to obtain visually acceptable shapes. However, instead of manually choosing this parameter, we can use additional information, which is greatly facilitated by the following observation:

Corollary 1 The estimation of λ can be performed after estimating µ and ν.

Proof This is a direct consequence of Prop. 3.

⊓ ⊔

This corollary justifies that µ and ν can be estimated before fixing ǫ: we can replace (33) by the following pair of problems, which must be solved sequentially:

         ( µ, ν) = argmin µ,ν TV(G(µ, ν, 1) ⊤ M 0 )    λ = argmin λ TV(G( µ, ν, λ) ⊤ M 0 ) s.t. λ ǫ (39)
where ǫ can then be chosen after estimating µ and ν. It can be stated that:

Lemma 1 The estimation of λ corresponds to the resolution of the classical bas-relief ambiguity.

Proof As we have

G( µ, ν, λ) ⊤ M 0 = G(0, 0, λ) ⊤ (G( µ, ν, 1) ⊤ M 0 ) (40)
we are now dealing with GBR matrices G(0, 0, λ) which will transform the field M 1 = G( µ, ν, 1) ⊤ M 0 : this particular GBR ambiguity with µ = ν = 0 corresponds to a bas-relief ambiguity on the field M 1 [START_REF] Belhumeur | The Bas-Relief Ambiguity[END_REF].

⊓ ⊔

The influence of λ is clear: λ → ∞ will "flatten" the depth function u while increasing the values of the albedo ρ, whereas λ → 0 will give "peaky" shapes with lower albedo. Indeed, according to [START_REF] Chambolle | A first-order primaldual algorithm for convex problems with applications to imaging[END_REF], the transformed depth ū writes ū(x, y) = u(x, y) λ

and, as shown in [START_REF] Belhumeur | The Bas-Relief Ambiguity[END_REF], the transformed albedo writes ρ(x, y) = ρ(x, y)

λ 2 + ∇u 2 1 + ∇u 2 1/2 (42)
Further discussion on the link between the geometry of the object and its albedo under bas-relief transformations can be found in [START_REF] Koenderink | The generic bilinear calibration-estimation problem[END_REF]. We show in Fig. 9 the effects of such a transformation on the shape and albedo of Beethoven's bust, estimated from the images of Fig. 1. As stated in Section 3, one possible solution involving total variation to solve for this last parameter is to minimize the total variation of the albedo:

(a) (b) (c) (d) (e) (f)
λ = argmin λ TV( G(0, 0, λ) ⊤ M 1 ) (43)
This approach is very similar to that of Alldrin et al. [START_REF] Alldrin | Resolving the Generalized Bas-relief Ambiguity by Entropy Minimization[END_REF] who minimize the entropy of the albedo distribution, nevertheless here only λ has to be estimated. However, we already stated in Section 3 that in general this approach gives poor results, as real-world albedos are usually not sufficiently regular, and moreover the cumbersome expression of the total variation of the albedo leads to numerical difficulties. Thus, it is better to use another method for estimating λ. In the experiments, we solve this residual ambiguity by adapting the constant light magnitude constraint [START_REF] Hayakawa | Photometric stereo under a light-source with arbitrary motion[END_REF]49] to the case of the bas-relief ambiguity, i.e. we assume that all the light vectors S i , i ∈ [1, m], have the same norm S 0 . When in addition µ and ν are unknown, it is shown in [49] that this constraint requires solving a constrained least-squares problem which is not trivial.

Here, µ and ν being already estimated, the problem is much easier. Let us S 1 = G( µ, 1) -1 S 0 . The constant light magnitude constraint reads:

∀i ∈ [1, m], (G(0, 0, λ) -1 S 1,i ) ⊤ (G(0, 0, λ) -1 S 1,i ) = S 2 0 (44)
where S 1,i refers to the i-th column of S 1 , S 0 is the unknown magnitude and λ is the bas-relief parameter. This is nothing more than a linear system with two unknowns (1/λ 2 , S 2 0 ) and m equations, which can be solved using Moore-Penrose pseudo-inverse. We can thus solve for λ 2 using this method, then deduce λ, since we assume λ > 0 for convexity. Note that this system might not admit a real solution (as the estimation of λ 2 might be negative), which is not stated in [START_REF] Hayakawa | Photometric stereo under a light-source with arbitrary motion[END_REF]49]. Such problems might occur for example if the distribution of the light magnitudes is far away from a uniform distribution: in such cases, other constraints should be considered. It is noticeable that the proposed method for estimating µ and ν is general, and can be used to simplify any method of resolution of the GBR. In Section 5, we adapt the Minimum-Entropy method [START_REF] Alldrin | Resolving the Generalized Bas-relief Ambiguity by Entropy Minimization[END_REF] when the constant light magnitude constraint is not valid: the discrete search on all the possible values of µ, ν and λ proposed in [START_REF] Alldrin | Resolving the Generalized Bas-relief Ambiguity by Entropy Minimization[END_REF] is reduced to a discrete search on the possible values of λ, which considerably accelerates the method.

An Alternative Approach: Depth Regularization

As will be shown in Section 5, the estimation of µ and ν proposed in Section 4.2, coupled with an efficient resolution of the bas-relief ambiguity, reaches state-ofthe-art results, and is considerably faster compared to other methods. It has however, in practice, the drawback of involving numerical approximations of the partial derivatives of M by finite differences. We can avoid computing these derivatives by regularizing the depth map u. Indeed, let us state that:

1. Formulation [START_REF] Mecca | Perspective Photometric Stereo with Shadows[END_REF] requires the knowledge of J(M 0 ). This involves necessary numerical approximations of ∇M 0 x , ∇M 0 y and ∇M 0 z in the implementation, which might occasionally be a source of problems, depending on the presence of discontinuities in the depth function. As such discontinuities are present in real-world objects, it might be necessary to apply Gaussian filtering on M 0 so as not to bias the results (the same "trick" is used during the integrability enforcement, as advised in [START_REF] Alldrin | Resolving the Generalized Bas-relief Ambiguity by Entropy Minimization[END_REF]). 2. We showed in Section 3 that the GBR has an impact on the variations of the albedo, the depth and its first order derivatives, however we did not explicitely use information about the depth.

Thus, we propose an alternative approach involving total variation of the depth u, which, as we will show, does not involve numerical approximations of ∇u by finite differences, as ∇u depends only on the known quantities p 0 = -M 0 x /M 0 z , q 0 = -M 0 y /M 0 z and on the GBR parameters.

Let us recall that the depth function u is assumed to be piecewise twice differentiable, so the integrability condition ∂ 2 u/∂x∂y = ∂ 2 u/∂y∂x is implicitely satisfied almost everywhere, and problem (31) can be reformulated in terms of depth as the following problem:

   u = argmin u TV(u) s.t. ρ (-∂xu,-∂yu,1) √ ∇u 2 +1 S -I 2 = min ( 45
)
Since this is as complicated to directly solve as [START_REF] Koenderink | The generic bilinear calibration-estimation problem[END_REF], we adopt the same parametric formulation. Let us denote by u 0 (x, y) the depth map obtained by integrating a normal field N 0 estimated by the method of Yuille and Snow [49] (we do not need to explicitely perform this integration), by U 0 = (x, y, u 0 (x, y)) ⊤ the corresponding Monge patch and by g(µ, ν, λ) = 1 λ (-µ, -ν, 1) the vector depending on the GBR parameters, which corresponds to the third line of G(µ, ν, λ) -1 . Similarly to Prop. 2, we have:

Proposition 4 A solution of (45) is given by u = g( µ, ν, λ)U 0 ( 46
)
where

   ( µ, ν, λ) = argmin µ,ν,λ TV(g(µ, ν, λ)U 0 ) s.t. λ > 0 ( 47
)
Proof Simply remark that, according to [START_REF] Chambolle | A first-order primaldual algorithm for convex problems with applications to imaging[END_REF], a GBRmodified depth function reads:

ū = g(µ, ν, λ)U 0 ( 48
)
The rest of the proof is similar to that of Prop. 2. ⊓ ⊔

We consider u 0 and U 0 are implicitely obtained by integrating the normal field estimated by using [49]. We focus on:

( µ, ν, λ) = argmin µ,ν,λ TV(g(µ, ν, λ)U 0 ) (49)
As the gradient of the modified depth function ū = g(µ, ν, λ)U 0 is equal to:

∇ū = - 1 λ (M 0 x /M 0 z + µ, M 0 y /M 0 z + ν) ⊤ = - 1 λ (M 0 x + µM 0 z )/M 0 z (M 0 y + νM 0 z )/M 0 z ( 50
)
the estimation of λ is problematic, as it to λ → ∞. In practice, similarly to the case, we have to impose a constraint λ ǫ ′ , ǫ ′ > 0 to avoid this degenerated solution. As previously, ǫ ′ can be manually tuned so as to obtain satisfactory 3D-reconstructions, or λ can be estimated by using another constraint, as already discussed in Section 4.3, which is justified by the following lemma:

Lemma 2 The estimations of µ an ν do not depend on that of λ.

Proof Eq. ( 50) shows that

∇(g(µ, ν, λ)U 0 ) = 1 λ ∇(g(µ, ν, 1)U 0 ) (51)
and thus:

argmin µ,ν TV(g(µ, ν, λ)U 0 ) = argmin µ,ν TV(g(µ, ν, λ ′ )U 0 ) (52)
for any λ ′ > 0. Note also that Eq. ( 51) shows that TV(u) is inversely proportional to λ (cf. Fig. 7-b). ⊓ ⊔ This lemma justifies the choice of separating the estimations of (µ, ν) and that of λ: we can first fix λ = 1 (for example) and estimate µ and ν. Once again, this reduces the generalized bas-relief ambiguity to the basrelief ambiguity.

Thus, we estimate µ and ν by solving:

( µ, ν) = argmin µ,ν TV(g(µ, ν, 1)U 0 ) ( 53
)
The solution of (53) can be obtained once again using quasi-Newton methods, the derivatives of E ′ (µ, ν) = TV(g(µ, ν, 1)U 0 ) = TV(-µx -νy + u 0 (x, y)) reading:

             ∂ µ E ′ = Ω M 0 x + µM 0 z (M 0 x + µM 0 z ) 2 + (M 0 y + νM 0 z ) 2 dx ∂ ν E ′ = Ω M 0 y + νM 0 z (M 0 x + µM 0 z ) 2 + (M 0 y + νM 0 z ) 2 dx ( 54
)
To accelerate the estimation, we choose to start from the solution ( µ 0 , ν 0 ) of the Tikhonov regularization:

( µ 0 , ν 0 ) = argmin µ,ν Ω ∇(g(µ, ν, 1)U 0 ) 2 dx (55)
which easily gives

µ 0 = p0 ν 0 = q0 ( 56
)
where p0 and q0 stand for the means of p 0 = -M 0 x /M 0 z and of q 0 = -M 0 y /M 0 z . In fact, (55) results in forcing the mean value of N to be equal to [0, 0, 1] ⊤ . This simply signifies that, on average, the normals should be oriented toward the camera, which is quite reasonable. This has also been used as a prior in an Expectation-Maximization algorithm for calibrated photometric stereo in presence of outliers [START_REF] Verbiest | Photometric stereo with coherent outlier handling and confidence estimation[END_REF].

Moreover, note that in (54), only known terms M 0 x , M 0 y and M 0 z are involved: no numerical approximation of J(M 0 ) is needed.

GBR-consistency

It was recently emphasized in [START_REF] Papadhimitri | A closed-form, consistent and robust solution to uncalibrated photometric stereo via local diffuse reflectance maxima[END_REF] that an appropriate method for estimating the GBR parameters should be GBR-consistent, i.e. it should return the same estimates M and S, whatever the initial solution (M 0 , S 0 ) of ( 8). Papadhimitri and Favaro prove in [START_REF] Papadhimitri | A closed-form, consistent and robust solution to uncalibrated photometric stereo via local diffuse reflectance maxima[END_REF] that this property holds for the "Diffuse Maxima" method (DM) initially described in [START_REF] Favaro | A Closed-Form Solution to Uncalibrated Photometric Stereo via Diffuse Maxima[END_REF]. However, this is not the case for the "Minimum Entropy" method [START_REF] Alldrin | Resolving the Generalized Bas-relief Ambiguity by Entropy Minimization[END_REF], neither for the "Self calibrating" method [START_REF] Shi | Self-calibrating Photometric Stereo[END_REF]. The estimations of µ and ν previously discussed actually hold this property, and we can state the following result: Proof We present the proof for the TV-regularization of M , the other proof is similar.

Let (M 01 , S 01 ) and (M 02 , S 02 ) be two different solutions of [START_REF] Bresson | Fast Dual Minimization of the Vectorial Total Variation Norm and Applications to Color Image Processing[END_REF]. The problem of estimating a field M with minimal TV, given the initial field M 01 , is synonymous to solving:

   ( µ, ν, λ) = argmin µ,ν,λ TV(G(µ, ν, λ) ⊤ M 01 ) s.t. λ > 0 (57)
On the other hand, there are necessarily three parameters (µ 01 , ν 01 , λ 01 ) ∈ R 2 × R * + such that the equality M 01 = G(µ 01 , ν 01 , λ 01 ) ⊤ M 02 holds. From the form of a GBR, we have G(µ, ν, λ) ⊤ G(µ 01 , ν 01 , λ 01 ) ⊤ = G(µ 01 + µλ 01 , ν 01 + νλ 01 , λλ 01 ) ⊤ . Writing µ ′ = µ 01 + µλ 01 , ν ′ = ν 01 + νλ 01 and λ ′ = λλ 01 , the problem becomes:

   ( µ, ν, λ) = argmin µ ′ ,ν ′ ,λ ′ TV(G(µ ′ , ν ′ , λ ′ ) ⊤ M 02 ) s.t. λ ′ > 0 (58)
which is exactly the problem of estimating a field M with minimal TV, given an initial field M 02 . Thus, the solution in M is the same, whether we initialize the algorithm with M 01 or with M 02 .

⊓ ⊔

In this section we evaluate the accuracy of the two posed methods. Some authors [19, evaluate this accuracy by calculating the mean angular error between the estimated normal field and that estimated by calibrated photometric stereo (Eq. ( 4)), although it is advocated in [START_REF] Alldrin | Resolving the Generalized Bas-relief Ambiguity by Entropy Minimization[END_REF] that uncalibrated photometric stereo might give overall "better" results than calibrated photometric stereo, as the use of prior information might compensate for some inaccuracy in the light parameters. Thus, we evaluate our methods firstly on synthetic datasets, for which ground truth normals are known.

Results on Synthetic Images

To evaluate the accuracy of the reconstructions obtained by our methods, referred to as TV-M and TV-u, we use synthetic images of a vase. The images were created in the following way:

• The depth equation can be written as the equation of a polynom (this equation can be found for example in [START_REF] Durou | Numerical methods for shape-from-shading: A new survey with benchmarks[END_REF]), thus the analytical expression of the normals is known. The depth function u is shown in Fig. 10-a, and its derivatives p and q in Figs. We compared our reconstructions to those of [START_REF] Alldrin | Resolving the Generalized Bas-relief Ambiguity by Entropy Minimization[END_REF] and [START_REF] Favaro | A Closed-Form Solution to Uncalibrated Photometric Stereo via Diffuse Maxima[END_REF], for which Matlab codes are freely available on the authors' websites. All our codes were developed using Matlab too so as to provide a fair comparison. CPU times were measured on a personal I7 processor at 2.9GHz. Those methods solve the GBR in this way:

• The Minimum-Entropy method (ME) [START_REF] Alldrin | Resolving the Generalized Bas-relief Ambiguity by Entropy Minimization[END_REF] estimates the GBR by minimizing the entropy of the albedo. • The Diffuse-Maxima method (DM) [START_REF] Favaro | A Closed-Form Solution to Uncalibrated Photometric Stereo via Diffuse Maxima[END_REF] considers the brightest pixels as locally oriented towards the light source and thus deduces the GBR.

As these methods have different features, we are studying three cases here, depending on the way we create the albedo:

• With a uniform albedo (Fig. 11-a), the albedo has a very low entropy, and the maxima of intensity give the orientation of the light, so both [START_REF] Alldrin | Resolving the Generalized Bas-relief Ambiguity by Entropy Minimization[END_REF] and [START_REF] Favaro | A Closed-Form Solution to Uncalibrated Photometric Stereo via Diffuse Maxima[END_REF] should work perfectly. • With a piecewise uniform albedo (Fig. 11-b), the albedo still has a low entropy, but the maxima of intensity might correspond to high albedo values, so this case favours [START_REF] Alldrin | Resolving the Generalized Bas-relief Ambiguity by Entropy Minimization[END_REF]. • With a radial gradient albedo (Fig. 11-c), the entropy is important, and maxima of intensity might again correspond to high albedo values.

One image of each dataset is shown in the second line of Fig. 11. To evaluate the accuracy of the different methods, we can calculate the relative errors between the estimated GBR parameters and the ground truth parameters, as well as the angular errors between the estimated normals and the ground truth normals. All these results are reported in Table 1. and the proposed methods TV-M (without or with applying Gaussian filtering (G.f) before calculating J(M 0 )) and TV-u. We show the relative errors on the GBR parameters µ, ν and λ, as well as the mean angular error (MAE), the maximal angular error (Max-AE), the minimal angular error (Min-AE), the standard deviation of the angular error (σ AE ) on the normals, and the CPU time. Uniform refers to the synthetic vase with uniform albedo, P/w unif. to the vase with piecewise uniform albedo, and Rad. grad. to the vase with a radial gradient albedo. The proposed methods overcome state-of-the-art results, and are considerably faster. They can also deal with situations state-of-the-art methods usually fail on (smoothly varying albedo).

As expected, the ME and DM methods give satisfactory results in the first case, but fail if the albedo varies smoothly. On the contrary, our methods produce comparable results in all three cases, as they do not rely on a hypothesis on the albedo, nor on pixel identification. They are thus more general methods, at least for estimating µ and ν.

The ME method uses a global optimization to minimize the (non-convex) entropy of the albedo distribution, resulting in very high CPU times. The CPU time for the DM method is also totally unpredictable, since it depends on the number of local maxima found in the images. On the contrary, the TV methods solely involve simple convex optimization tools and are thus faster.

Finally, let us explain the differences between the results of TV-M and TV-u. As mentioned, TV-M involves numerical approximations, which is totally biased by the non-integrability of the vase on the boundaries. It is for that reason that TV-u performs better. It is possible to apply Gaussian filtering before calculating the numerical derivatives in [START_REF] Mecca | Perspective Photometric Stereo with Shadows[END_REF] and [START_REF] Onn | Integrability Disambiguates Surface Recovery in Two-Image Photometric Stereo[END_REF], which reduces the impact of these approximations and improves the results of the TV-M method. This modified method is referred to as TV-M (G.f).

Evaluation on Real-World Datasets

Datasets

To evaluate the accuracy of the reconstructions on realworld datasets, since no ground truth is available, we can only compare the results to those of calibrated photometric stereo. It should be recalled that this comparison is purely informative, since estimating the normal field by calibrated photometric stereo is subject to inaccuracy in the given light parameters. We used:

• Ten datasets of almost Lambertian objects. The Buddha, Cat, Horse, Owl and Rock datasets, courtesy of Dan Goldman and Steven Seitz4 , are composed of twelve 320 × 300 images (except Horse, which is composed of twelve 270×350 images). The Octopus, Redfish and Korean Doll datasets can be found on Neil Alldrin's homepage5 : they are composed of five images, respectively of size 321 × 281, 351 × 301 and 1321 × 521. The Doll dataset was presented in [START_REF] Favaro | A Closed-Form Solution to Uncalibrated Photometric Stereo via Diffuse Maxima[END_REF], and is composed of fifteen 405×250 images. Finally, Beethoven dataset is composed of three 256 × 256 images.

• The first twenty sets of images the Extended Yale DataFace B [START_REF] Georghiades | From Few to Many: Illumination Cone Models for Face Recognition under Variable Lighting and Pose[END_REF][START_REF] Lee | Acquiring Linear Subspaces for Face Recognition under Variable Lighting[END_REF]. These datasets originally contain 64 different illumination conditions: we selected only the 27 with lowest elevation angles, so as to discard the images the most corrupted by shadows. Each image of each dataset is 192 × 168.

Because real-world images in reality contain outliers such as shadows we apply the preprocessing from [START_REF] Wu | Robust Photometric Stereo via Low-Rank Matrix Completion and Recovery[END_REF], which involves a hyper-parameter C that we fixed, as advised in [START_REF] Favaro | A Closed-Form Solution to Uncalibrated Photometric Stereo via Diffuse Maxima[END_REF], to for the datasets with at least 12 images, and to 3 for the other datasets.

One preprocessed image of each dataset is shown in Fig. 12. However such a preprocessing cannot remove all the outliers: specular highlights are still visible in, for instance, the Owl dataset (see the eyes of the owl in Fig. 12).

Results

The relative errors on the GBR parameters, the mean angular errors on the normal fields and the CPU times can be found in Tables 2 and3. As can be seen, the proposed methods reach, in general, state-of-the-art results on those datasets, and considerably overcome the two others in terms of computation times.

For some datasets, no error is given: this corresponds to the cases where the resolution of (44) results in negative values of λ 2 . The constraint of uniform light magnitude is not adapted to such datasets. However, as the estimated values of µ and ν do not rely on this constraint, the estimation of λ can be achieved for example by adapting another method of resolution of the GBR to the resolution of the bas-relief ambiguity. Thus, we also report the results of estimating λ by ME, with µ and ν estimated by TV-u (this method is referred to as TV-u / ME). Unifying models [START_REF] Chandraker | Reflections on the Generalized Bas-Relief Ambiguity[END_REF] and (53), this method can be written:

     ( µ, ν) = argmin µ,ν TV(g(µ, ν, 1)U 0 ) λ = argmin λ E( G( µ, ν, λ) ⊤ M 0 ) ( 59
)
The discrete search for the minimization of the entropy is thus performed only with respect to λ, which results in considerably shorter computation times, compared to ME. The example of the Korean Doll dataset (cf. Table 2) is striking: the computation time is reduced from more than two minutes to a couple of seconds using TV-u / ME, while obtaining a very acceptable level of accuracy. This confirms that the proposed method for the estimation of µ and ν can be used to accelerate any method of GBR resolution.

Moreover, these tests empirically prove the validity of the "Minimum-TV" constraint for a wide variety of objects, since the obtained errors are comparable to state-of-the-art methods. However, as the proposed methods do not rely either on pixel identification or on minimizing non-convex functions, they are faster than other methods.

Finally, we show in Fig. 13 a side-view of a relighting of each 3D-model obtained after integrating the normals estimated with the TV method which gives the lowest MAE and warping the estimated albedo on the surface. We use three light sources (one in front of the object, one on the left side and one on the right side) and the Phong model [START_REF] Phong | Illumination for computer generated pictures[END_REF] for rendering.

Guidelines

Finally, we provide some guidelines to help the user choose the right method. To use the proposed methods, two choices have to be made:

• TV-regularization of M or TV-regularization of u:

this is linked to the estimation of µ and ν. • How to solve the residual bas-relief ambiguity, which is linked to the estimation of λ.

The TV-regularization of u should, in the most general case, give better results, as it does not rely on numerical differentiations of the initial field M 0 , which are biased in the presence of depth discontinuities (e.g., the Owl dataset). However, if both the surface and the albedo are smooth enough (e.g., the Horse and the B13 datasets), the TV-regularization of M will offer comparable or better results. However, it is hard to a priori ensure the absence of depth discontinuities.

Thus, to choose the appropriate method, we advise considering the following strategy:

• If no prior information is available, the best choice is probably to firstly try the TV-u method: the estimation of µ and ν does not involve any numerical approximation, and the estimation of λ using the constant light magnitude constraint is quite fast. • If the surface and the albedo seem smooth enough, the TV-M method may give overall better estimations of µ and ν. • Finally, if the constant light magnitude constraint is not satisfied, in both cases (TV-u or TV-M ) negative values of λ 2 will be obtained, or the surface will look either too flat or too peaky. In that case, other methods of estimation of λ should be considered. If the albedo is regular enough, one can easily adapt the Min. Entropy method to the estimation of λ only (we provide such an implementation in Matlab): this corresponds to the TV-u / ME method.

Other methods [START_REF] Favaro | A Closed-Form Solution to Uncalibrated Photometric Stereo via Diffuse Maxima[END_REF][START_REF] Shi | Self-calibrating Photometric Stereo[END_REF] would also be easy to adapt. 1.

Conclusion and Future Work

Contributions. In this paper, we showed that the properties of total variation were highly useful in the resolution of the uncalibrated photometric stereo problem, since estimating the parameters of the GBR transformation as the minimizers of the total variation of either the M field or the depth function u reduces the GBR ambiguity to the classical bas-relief ambiguity, which is considerably easier to solve. We experimentally proved that the proposed methods reached state-of-theart results in terms of accuracy, with much shorter computation times, thanks to the parametric formulation of the original problem [START_REF] Koenderink | The generic bilinear calibration-estimation problem[END_REF], which we converted into the problem of minimizing a convex function depending only on 3 parameters.

Future prospects. Despite the fact that very satisfactory results were obtained using the proposed methods, their justification remains empirical: a statistical study of the distributions of ρ, u and ∇u over a wide variety of objects would probably help understanding why total variation offers such satisfactory results. All the work presented here assumes the light sources considered are directional, so as to obtain a linear formulation of the problem. In practice, directional lighting is difficult to ensure, and point-light sources should be considered. The problem becoming nonlinear, the GBR ambiguity would disappear. However, the problem would become numerically much more difficult, so efficient strategies should be considered.

Finally, for a robust estimation, it would be necessary to consider outliers such as noise, specularities and shadows, as they induce a bias in the estimation. Indeed, as the preprocessing [START_REF] Wu | Robust Photometric Stereo via Low-Rank Matrix Completion and Recovery[END_REF] is purely driven by the rank of the matrix, it does not take spatial considerations into account and thus gives poor results in presence of noise or strong outliers, therefore some artifacts in the reconstructions are still visible (see the eye of the owl in Fig. 13). However, such artifacts could easily be removed by "smoothing" the normal field: the model (30) could be used as a post-processing method after the estimation of M and before the integration step. To ensure that the normal fields are "almost everywhere" integrable, we could switch the L 2 minimization curl M 2 = min in [START_REF] Horovitz | Depth from gradient fields and control points: Bias correction in photometric stereo[END_REF] to a L 1 minimization curl M 1 = min (this could be efficiently achieved by using iteratively reweighted least-squares). Proceeding so would actually result in an improvement of the method of Reddy et al. [START_REF] Reddy | Enforcing integrability by error correction using L1minimization[END_REF], who are concerned in enforcing integrability on the normal field, but not in the particular case of photometric stereo, for which not only the normal field, but also the albedo has to be corrected.

Fig. 1 :

 1 Fig. 1: m = 3 images, of size 256 × 256, of a Beethoven's bust, taken under 3 different illuminations 2 .

Fig. 2 :

 2 Fig. 2: Results of applying photometric stereo to the three images of Fig. 1: (a) ρ; (b) p; (c) q. (d) Relighting of the reconstructed surface. White pixels in (a)-(b)-(c) indicate high values.

Fig. 3 :

 3 Fig. 3: Effects of a GBR of parameters µ = 1/2, ν = 1 and λ = 1/3 on: (a) the albedo, (b)-(c) the normals and (d) the shape of Fig. 2. The scales are the same as in Fig.2, so as to illustrate the drifts in all the three directions. It is to be noted that the albedo (a) is spread, compared to that of Fig.2-a: this led to the Minimum-Entropy algorithm[START_REF] Alldrin | Resolving the Generalized Bas-relief Ambiguity by Entropy Minimization[END_REF].

Fig. 4 :

 4 Fig. 4: Three different albedo configurations and the corresponding histograms: (a) has a significant entropy value, cf. histogram (d); (b) and (c) have the same slight entropy, cf. histograms (e) and (f), however the spatial distributions of the albedo are different. Total variation would tend to favour distributions such as (b).

Fig. 5 :

 5 Fig. 5: The effects of a GBR transformation on the variations of the albedo, of the depth and of its derivatives. (a) One of the input images. (b) The y = 100 straight line. (c) Ground truth 3D-model. (d) GBR-transformed 3D-model.We plotted the values of the ground truth and the GBR-transformed values of (e) ρ, (f) u, (g) p and (h) q along the y = 100 straight line, with µ = 1/2, ν = 1 and λ = 1/3. Obviously, the GBR-transformed estimations have higher total variations.

Fig. 6 :

 6 Fig. 6: The effects of µ and ν on the total variations of (a) the albedo ρ, (b) the depth function u and (c) the M field. The three total variations seem to be minimal for values of µ and ν close to the ground truth values.

Fig. 7 :

 7 Fig. 7: The effects of λ on the total variations of (a) the albedo ρ, (b) the depth function u and (c) the M field. Only the total variation of the albedo has a minimum inside ]0; +∞[, whereas TV(u) is inversely proportional to λ and TV(M ) is proportional to λ. Both these last observations are proven in Section 4.

Fig. 8 :

 8 Fig. 8: The total variation of the albedo ρ, as a function of λ, on the Doll dataset. The value of λ minimizing TV(ρ) is quite far away from the ground truth value 1.

Fig. 9 :

 9 Fig. 9: Effects of a bas-relief transformation. First line: reconstructed shapes with (a) λ = 1/2, (b) λ = 1 (original shape), and (c) λ = 2. Second line: albedo with (d) λ = 1/2, (e) λ = 1, and (f) λ = 2.

Proposition 5

 5 Both the TV methods described in Sections 4.2 and 4.4 are GBR-consistent.

  10-b and 10-c. • We created a 3 × m light matrix S, choosing m = 22 light vectors having the same intensity (so as to solve the residual bas-relief ambiguity). • We manually created an albedo map, thus obtaining a ground truth M field (see below). • Eventually, the 256 × 256 images were generated using the Lambertian model: I = max(M ⊤ S, 0) (the max operator models the self-shadows).

Fig. 10 :

 10 Fig. 10: Ground truth functions (a) u, (b) p and (c) q for the synthetic vase dataset.

Fig. 11 :

 11 Fig. 11: First line: ground truth albedo maps: (a) uniform; (b) piecewise uniform; (c) radial gradient. Second line: one image of each dataset, corresponding to each different albedo (the same lighting was used to generate these three images).

Fig. 12 :Fig. 13 :

 1213 Fig. 12: One preprocessed image of each real-world dataset used in the experiments.

Table 1 :

 1 Comparison of the results of the Minimum-Entropy method (ME), the Diffuse-Maxima method (DM)

		|µ-µ| |µ|	|ν-ν| |ν|	|λ-λ| |λ|		Max-AE (deg.) Min-AE (deg.) σ AE (deg.) CPU (s)
	ME	1.54	0.00	0.05	0.63	16.41	0.03	1.10	26.41
	DM	5.36	0.28	0.01	0.98	16.04	0.01	0.95	2.79
	TV-M	93.72	0.27	0.57	10.09	24.77	0.03	3.14	0.13
	TV-M (G.f)	0.61	0.02	0.07	0.57	16.75	0.00	1.28	0.26
	TV-u	0.61	0.02	0.07	0.57	16.75	0.00	1.28	0.57
	P/w unif.
	ME	1.19	0.02	0.04	0.59	16.35	0.00	1.09	28.54
	DM	24.11	1.59	0.11	4.82	14.97	0.06	1.25	3.34
	TV-M	53.38	0.11	0.36	7.42	21.53	0.11	2.11	0.07
	TV-M (G.f)	0.50	0.02	0.07	0.56	16.77	0.00	1.27	0.26
	TV-u	0.50	0.02	0.07	0.56	16.77	0.00	1.27	0.50
	Rad. grad.
	ME	8.01	18.22	0.21	15.59	18.89	2.03	3.64	30.82
	DM	8.77	20.74	0.10	15.53	18.40	3.03	3.45	0.56
	TV-M	102.61	0.89	0.20	4.43	18.76	0.07	1.35	0.06
	TV-M (G.f)	2.70	0.03	0.07	0.75	16.53	0.01	1.26	0.15
	TV-u	2.41	0.02	0.07	0.75	16.53	0.01	1.26	0.23

Table 2 :

 2 Comparison of the results of the ME, DM and TV methods, on the first twenty datasets of the extended Yale DataFace B. Note that B14 does not exist in this database. Units are the same as in Table1.

		|µ-µ| |µ|	|ν-ν| |ν|	|λ-λ| |λ|	MAE	CPU
	Buddha
	ME	0.81	0.22	0.06	6.38	13.87
	DM	0.07	0.06	0.12	4.99	1.91
	TV-M	0.91	0.35	0.21	10.92	0.06
	TV-u	1.03	0.51	1.34	18.85	0.03
	TV-u / ME	1.03	0.51	0.00	14.79	0.35
	Cat
	ME	0.48	0.49	0.17	14.37	17.82
	DM	0.95	0.05	0.08	5.37	0.57
	TV-M	1.00	0.10	0.27	6.16	0.06
	TV-u	0.13	0.54	5.94	35.57	0.04
	TV-u / ME	0.13	0.54	0.18	16.04	0.45
	Horse
	ME	11.33	0.28	0.09	10.88	14.21
	DM	14.25	0.20	0.03	4.80	1.25
	TV-M	11.64	0.11	0.21	6.68	0.07
	TV-u	34.83	0.34	0.64	18.25	0.03
	TV-u / ME	34.83	0.34	0.13	15.59	0.34
	Owl
	ME	0.41	0.63	0.25	15.67	22.71
	DM	0.35	0.00	0.07	6.63	0.43
	TV-M	0.50	0.23	0.13	8.42	0.08
	TV-u	0.15	0.13	0.08	5.54	0.05
	TV-u / ME	0.15	0.13	0.15	7.54	0.54
	Rock
	ME	14.77	0.75	0.23	18.37	32.30
	DM	5.35	0.39	0.05	11.61	0.61
	TV-M	3.13	0.88	-	-	0.08
	TV-u	9.43	1.19	-	-	0.07
	TV-u / ME	9.43	1.19	0.19	27.56	0.76
	Octopus
	ME	0.06	0.07	0.07	3.06	22.49
	DM	0.10	0.16	0.19	6.64	0.26
	TV-M	0.03	0.14	-	-	0.07
	TV-u	0.11	0.29	-	-	0.07
	TV-u / ME	0.11	0.29	0.28	9.41	0.67
	Redfish
	ME	7.06	0.06	0.03	4.38	21.52
	DM	8.85	0.18	0.05	5.60	0.79
	TV-M	20.83	0.06	3.48	25.92	0.10
	TV-u	3.28	3.75	2.93	26.07	0.05
	TV-u / ME	3.28	3.75	0.17	14.17	0.61
	Korean Doll
	ME	0.04	0.19	0.08	3.39	151.99
	DM	0.20	0.80	0.27	17.74	2.28
	TV-M	0.11	0.46	-	-	0.49
	TV-u	0.09	0.50	-	-	0.34
	TV-u / ME	0.09	0.50	0.15	10.36	3.87
	Doll
	ME	0.05	2.64	0.47	26.24	17.13
	DM	0.03	0.81	0.34	12.15	0.59
	TV-M	0.05	1.81	-	-	0.06
	TV-u	0.04	1.82	-	-	0.03
	TV-u / ME	0.04	1.82	0.44	19.84	0.35
	Beethoven
	ME	10.50	6.69	0.57	63.83	16.39
	DM	0.02	0.25	0.41	10.01	1.31
	TV-M	0.15	0.34	0.65	16.44	0.05
	TV-u	0.06	0.31	0.47	13.70	0.03
	TV-u / ME	0.06	0.31	0.07	13.30	0.32

Table 3 :

 3 Comparison of the results of the ME, DM and TV methods, on ten almost Lambertian objects. Units are the same as in Table

The Matlab codes corresponding to the proposed methods are available online at: http://ubee.enseeiht.fr/photometricstereo/

However, as mentioned in the perspectives, the model[START_REF] Horovitz | Depth from gradient fields and control points: Bias correction in photometric stereo[END_REF] could be useful to post-process the results.

http://courses.cs.washington.edu/courses/cse455/10wi/projects/

http://vision.ucsd.edu/ ~nalldrin/research/