Yvain Quéau

Roberto Mecca

Jean-Denis Durou

Xavier Descombes

Photometric Stereo with Only Two Images: A Theoretical Study and Numerical Resolution

Keywords: 3D-reconstruction, shape-from-shading, photometric stereo, PDEs, numerical analysis, optimization, graph cut

scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In the computer vision field, 3D-shape reconstruction using digital images as input data has gained a growing importance. Interest in this task has increased even more since most mass digital devices have been equipped with cameras. Based on more than thirty years of research, such devices are $ Author version. Publisher version's DOI: 10.1016/j.imavis.2016.11.006

Email address: yvain.queau@tum.de (Yvain Quéau) 1 Marie Curie fellow of the Istituto Nazionale di Alta Matematica, Italy.

potentially convertible into 3D-scanners without any hardware correction. Among all the photographic 3D-reconstruction techniques, we focus in this work on shape-from-shading (SFS) and photometric stereo (PS), which exploit shading information when one (SFS) or several (PS) sources illuminate the observed object. For a comprehensive overview on these techniques, see the reference book [START_REF]Shape from Shading[END_REF] by Horn and Brooks, but also [START_REF] Durou | Numerical Methods for Shape-fromshading: A New Survey with Benchmarks[END_REF] and [START_REF] Vogiatzis | Practical 3D Reconstruction Based on Photometric Stereo[END_REF][START_REF] Ackermann | A survey of photometric stereo techniques[END_REF] for up-to-date surveys on SFS and PS, respectively.

Many articles enlighted the impossibility of avoiding any ambiguity while retrieving the shape from a single image, as in the SFS problem [5]. This impossibility arises from the difficulty met in distinguishing the concave from the convex surfaces. The most natural way to solve this problem is to use more than one image: Woodham showed in [6] that three is the minimum number of images to ensure well-posedness of the PS problem.

This work focuses on the intermediate case when only two images are taken into account (this specific situation will be referred to as PS2). Besides being particularly interesting for dedicated applications as single-day outdoor PS from sun light [START_REF] Sato | Reflectance analysis under solar illumination[END_REF][START_REF] Jung | One-day outdoor photometric stereo via skylight estimation[END_REF], the PS2 problem can be seen as the degenerative case of lack of information from the three-source PS problem due to shadows [START_REF] Hernández | Overcoming shadows in 3-source photometric stereo[END_REF]. In this view, we provide working tools aimed at solving the underlying ambiguities, derived from a theoretical study.

There exist a combinatorial number of normal fields which are solutions of the PS2 problem. Exhaustive search can be carried out among these normal fields, in order to find the one which best satisfies a smoothness constraint [START_REF] Onn | Integrability Disambiguates Surface Recovery in Two-Image Photometric Stereo[END_REF]. Alternatively, one may resort to the differential approach of PS, which implicitly enforces smoothness. A meaningful solution of the resulting PDE can be obtained either by specifying an explicit boundary condition [START_REF] Mecca | Unambiguous Photometric Stereo Using Two Images[END_REF], or by resorting to regularization [START_REF] Hernández | Overcoming shadows in 3-source photometric stereo[END_REF]. Unfortunately, knowledge of the surface on the boundary is rarely available, and regularization techniques come along with parameters to tune, which might be tedious.

We put forward a new method for solving the PS2 problem, which is based on the non-differential approach of PS (normal estimation). We assume that the setup consists of orthographic viewing geometry, parallel and uniform lighting, as well as Lambertian reflectance. Under such assumptions, and provided that the albedo is known, it is shown in this paper that the number of possible solutions can be predicted beforehand. It is then demonstrated that exhaustive search of the "best" normal field can be recast as a binary labelling problem, efficiently solvable by resorting to the graph cut algorithm. In contrast with existing methods, the proposed one requires neither knowing a boundary condition, nor tuning any parameter.

In the following, we first review the general equations of SFS, PS and PS2 in Section 2. In Section 3, we compare the differential and non-differential formulations of PS2. We show why the PS2 problem usually has a unique solution in Section 4. In Section 5, a practical graph cut-based [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF] algorithm to compute this solution is introduced, and it is applied in Section 6 to the problem of three-source PS with shadows.

2. From Shape-from-shading to Photometric Stereo

Shape-from-shading

To fully describe the problem we are interested in, let us first recall some features of the SFS problem. We attach to the camera a 3D-Cartesian coordinate system xyz, so that xy coincides with the image plane and z with the optical axis. Under the assumption of orthographic projection, the visible part of a surface is a graph z = u(x, y). It is well known that the SFS problem is modeled by the image irradiance equation [START_REF]Shape from Shading[END_REF]: R(n(x, y)) = I(x, y)

where I(x, y) is the graylevel at image point (x, y), and the reflectance function R(n(x, y)) gives the value of the light re-emitted by the surface as a function of its orientation i.e., of the unit-length outgoing normal n(x, y) to the surface at surface point [x, y, u(x, y)] . The unknown depth u has to be reconstructed on a compact domain Ω ⊂ R 2 called the reconstruction domain.

Let us consider a unique parallel and uniform light beam whose direction is indicated by the unit-length vector s = [s 1 , s 2 , s 3 ] = [s , s 3 ] ∈ R 3 , and whose intensity is denoted by ψ. Assuming the observed object has purely diffuse reflection, and ignoring shadows, Eq. ( 1) can be written as follows:

ρ(x, y) ψ s • n(x, y) = I(x, y) (2)
where ρ(x, y) ∈ [0, 1] is the albedo. In fact, this equality is nothing more than a relation of proportionality. Knowing that the vectors s and n(x, y) have unit-length, and assuming that ψ is a constant factor, it seems justified to rewrite Eq. (2) as a real equality:

ρ(x, y) s • n(x, y) = I(x, y) (3)
where I(x, y) ∈ [0, 1] should now be considered as the normalized graylevel. Eq. ( 3) is a particular non-differential formulation (amongst many others) of the SFS problem. Once the normal field n has been estimated, it has to be integrated. This means that the following equation in u has to be solved [START_REF] Durou | Integration of a Normal Field in the Presence of Discontinuities[END_REF]:

n(x, y) = 1 1 + ∇u(x, y) 2 -∇u(x, y) , 1 (4)
where ∇u(x, y) = [∂ x u(x, y), ∂ y u(x, y)] denotes the gradient of u(x, y). From Eqs. ( 3) and ( 4), we get the following differential formulation of SFS:

ρ(x, y) -s • ∇u(x, y) + s 3 1 + ∇u(x, y) 2 = I(x, y) (5)
which is a first-order nonlinear PDE of the Hamilton-Jacobi type. We refer the interested reader to the survey presented in [START_REF] Durou | Numerical Methods for Shape-fromshading: A New Survey with Benchmarks[END_REF] for a presentation of recent results on the eikonal equation, which follows from (5) when s = [0, 0, 1] .

Photometric Stereo

Even if s is known, SFS is ill-posed without any additional knowledge on the surface to be reconstructed. In most papers on SFS, the albedo ρ(x, y) is supposed to be known, but this is still not enough to make the problem well-posed. The simplest way to overcome SFS ill-posedness is to use m ≥ 2 images taken from the same point of view, illuminated by m light sources (s i , ψ i ), i ∈ [1, m]. This new problem is called photometric stereo (PS). The classical resolution of PS is based on a local estimate of the outgoing unitlength normal to the surface [6]. For a Lambertian surface, the non-differential formulation of PS consists in solving a system of m equations of type (3):

ρ(x, y) s i • n(x, y) = I i (x, y), i ∈ [1, m] (6)
As for SFS, this formulation requires that Eq. ( 4) is solved afterwards. From a theoretical point of view, Eq. ( 4) admits a solution in u only if the estimated normal field is integrable [START_REF] Frankot | A Method for Enforcing Integrability in Shape from Shading Algorithms[END_REF] (cf. Section 4.5). Due to estimation errors, this is rarely the case in practice, hence projection of the estimated normal field on the space of integrable fields must be achieved, resorting for instance to Fourier analysis [START_REF] Frankot | A Method for Enforcing Integrability in Shape from Shading Algorithms[END_REF] or to variational methods [START_REF] Durou | Integration of a Normal Field in the Presence of Discontinuities[END_REF]. Alternatively, one can directly try to estimate the "most integrable" normal field. This is the approach that is followed in Section 5.

Of course, a differential formulation of PS also exists, which aims at solving a system of m nonlinear PDEs of type (5):

ρ(x, y) -s i • ∇u(x, y) + s i 3 1 + ∇u(x, y) 2 = I i (x, y), i ∈ [1, m] (7)
In the usual case, denoted PS3, where m ≥ 3 non-coplanar distant calibrated light sources are used [6], System (6) reduces to a full-rank linear system in m(x, y) = ρ(x, y) n(x, y) ∈ R 3 . Solving this system has several advantages, compared to SFS: it is well-posed and can be locally solved, thus parallelized. Furthermore, the albedo no longer has to be known.

We may wonder whether the differential formulation [START_REF] Sato | Reflectance analysis under solar illumination[END_REF] would really be pertinent for PS3. The main advantage of solving [START_REF] Sato | Reflectance analysis under solar illumination[END_REF] is that integrability is implicitly ensured, unlike solving (6), knowing that the lack of integrability of n complicates the resolution of Eq. ( 4) [START_REF] Frankot | A Method for Enforcing Integrability in Shape from Shading Algorithms[END_REF]. However, the problem [START_REF] Sato | Reflectance analysis under solar illumination[END_REF] has two drawbacks: it is nonlinear and cannot be solved locally [START_REF] Mecca | Uniqueness and Approximation of a Photometric Shape-from-Shading Model[END_REF].

In this paper, we focus on the resolution of PS when the linear system (6) is not full-rank. In such cases, the non-differential formulation does not have as many advantages as for PS3, and the differential formulation might be worthwhile, at least because it is better-posed since the integrability constraint is implicitly satisfied. Indeed, practical solutions to the rank-deficient PS problem use this differential formulation. Yet, as discussed in Section 3.2, differential approaches have to resort either to a boundary condition (which is rarely available) or to regularization (which requires parameter tuning). The solution presented in Section 5, which is based on the non-differential formulation, involves neither boundary condition nor parameter tuning.

Scope of our Work: Photometric Stereo Using Two Images

The scope of our work is the intermediate case m = 2 of PS between m = 1 (SFS) and m ≥ 3 (PS3), denoted PS2. Such a problem arises in real-time 3D-reconstruction of non-rigid objects [START_REF] Hernández | Nonrigid Photometric Stereo with Colored Lights[END_REF]. This application can be carried out using PS3, an RGB sensor and three colored light sources. Yet, when one of the light source is occluded, only two of the color channels provide meaningful shading clues, and hence a PS2 problem must be solved. The PS2 problem is also very similar to the case m ≥ 3 with coplanar light vectors s i . A concrete example to be mentioned is the case when a scene is illuminated by the sun [START_REF] Sato | Reflectance analysis under solar illumination[END_REF]. In both cases, the linear system (6) is no longer full-rank in m(x, y) = ρ(x, y) n(x, y).

The PS2 problem has been addressed in few papers. Onn and Bruckstein prove in [START_REF] Onn | Integrability Disambiguates Surface Recovery in Two-Image Photometric Stereo[END_REF] that the determination of the normal is a priori ambiguous, and demonstrate how to eliminate this ambiguity using the integrability constraint. Yang et al. study the problem in the particular case of convex objects [START_REF] Yang | Two image photometric stereo method[END_REF]. In [START_REF] Sato | Reflectance analysis under solar illumination[END_REF], Sato and Ikeuchi use the resolution method designed by Onn and Bruckstein to solve the PS problem using m ≥ 3 images under solar illumination, which partially brings us back to the PS2 problem [18]. In [19], Kozera reaches the same conclusions as Onn and Bruckstein, by making an analytical resolution of the differential formulation. From 1995 and for over ten years, only Ikeda addressed the PS2 problem [20]. However, he essentially considers the second image as a means to more accurately solve the SFS problem. More recently, the problem of outdoor PS has been reexplored in [START_REF] Jung | One-day outdoor photometric stereo via skylight estimation[END_REF]21,[START_REF] Ackermann | Photometric Stereo for Outdoor Webcams[END_REF]. Finally, let us quote a preliminary version [START_REF] Mecca | Unambiguous Photometric Stereo Using Two Images[END_REF] of our work.

One could wonder whether the PS2 problem is not purely formal, considering that the actual trend is to deal with far more than two images, as in [23] where a full video is used, or in [21,[START_REF] Ackermann | Photometric Stereo for Outdoor Webcams[END_REF] where timelapse images sequences, acquired over several years, are considered. Apart from the intrinsic interest we can find in studying the number of solutions to the PS2 problem, PS in the presence of shadows has been studied in several recent papers [START_REF] Hernández | Overcoming shadows in 3-source photometric stereo[END_REF]24,25,26,[START_REF] Verbiest | Photometric stereo with coherent outlier handling and confidence estimation[END_REF][START_REF] Wu | Robust Photometric Stereo via Low-Rank Matrix Completion and Recovery[END_REF][START_REF] Ikehata | Photometric Stereo Using Sparse Bayesian Regression for General Diffuse Surfaces[END_REF]. When m > 3, the graylevels which lie inside shadows are considered as outliers and left out of the estimation [24, 25,26] or dealt with in a robust estimation process such as the Expectation-Maximization algorithm [START_REF] Verbiest | Photometric stereo with coherent outlier handling and confidence estimation[END_REF]. Most recent techniques assume that such outliers are sparsely distributed in the images, and hence design sparsity-enhancing algorithms [START_REF] Wu | Robust Photometric Stereo via Low-Rank Matrix Completion and Recovery[END_REF][START_REF] Ikehata | Photometric Stereo Using Sparse Bayesian Regression for General Diffuse Surfaces[END_REF]. When m = 3, Hernández et al. show in [START_REF] Hernández | Overcoming shadows in 3-source photometric stereo[END_REF] that implementing a specific treatment for the twice-lit points i.e., for those where the PS2 problem may arise, improves the 3D-reconstruction accuracy. The application of our study will precisely be to improve the 3D-reconstruction in such shadow areas (see Section 6).

Photometric Stereo Using Two Images: A Theoretical Study

Non-differential Formulation of the PS2 Problem

Taking m = 2 in (6), the non-differential formulation of PS2 is written as:

     ρ(x, y) s 1 1 n 1 (x, y) + s 1 2 n 2 (x, y) + s 1 3 n 3 (x, y) = I 1 (x, y) ρ(x, y) s 2 1 n 1 (x, y) + s 2 2 n 2 (x, y) + s 2 3 n 3 (x, y) = I 2 (x, y) n 1 (x, y) 2 + n 2 (x, y) 2 + n 3 (x, y) 2 = 1 (8)
If ρ(x, y) is assumed to be known, the problem (8) consists at each point (x, y) in a nonlinear system of three equations with three unknowns (n 1 , n 2 , n 3 ). The nonlinearity of the third equation could give rise to a non-unique solution, as observed by Ikeuchi and Horn in [START_REF] Ikeuchi | Numerical Shape from Shading and Occluding Boundaries[END_REF]: "Naturally, the above nonlinear equations [...] may have more than one solution, in which case additional information (such as a third image) may be needed to find a unique answer". Indeed, most of the time, [START_REF] Jung | One-day outdoor photometric stereo via skylight estimation[END_REF] does not have a unique but two solutions.

To solve [START_REF] Jung | One-day outdoor photometric stereo via skylight estimation[END_REF], we only consider the set S of twice-lit unit-length normals n i.e., those where s 1 • n > 0 and s 2 • n > 0. This set is an open part of Gaussian sphere G limited by two planes containing the origin (see Fig. 1-a):

π 1 : s 1 1 n 1 + s 1 2 n 2 + s 1 3 n 3 = 0 π 2 : s 2 1 n 1 + s 2 2 n 2 + s 2 3 n 3 = 0 (9)
Planes π 1 and π 2 are orthogonal, respectively, to s 1 and s 2 . At each twicelit point (x, y) characterized by the graylevels I 1 (x, y) and I 2 (x, y), and by the supposedly known albedo ρ(x, y) = 0, System (8) admits two solutions, denoted n + (x, y) and n -(x, y), which are the intersections of S and of two planes which are obtained by translating π 1 and π 2 of I 1 (x, y)/ρ(x, y) and I 2 (x, y)/ρ(x, y) in the directions of s 1 and s 2 , respectively. See Fig. 1-b for a geometrical interpretation of n + (x, y) and n -(x, y). The analytical formulae for them will be derived in Section 5.1. Let us denote by π the plane supported by s 1 and s 2 and containing the origin of G. It is obvious that n + (x, y) and n -(x, y) are symmetric with respect to π, for any (x, y) ∈ Ω. If I 1 (x, y) and I 2 (x, y) exactly match the Lambertian model (3), the non-differential problem (8) admits either two solutions, or one if n + (x, y) = n -(x, y). But, if this model is not exactly satisfied, which may happen with real images, there may be some points (x, y) without any exact solution. The way to effectively handle these different cases will be specified in Section 5.1.

π 2 π 1 G s 2 s 1 S E x z y E s 2 z y x s 1 (a) (b)
Figure 1: (a) The set S of twice-lit normals, emphasized in red, is a part of Gaussian sphere G limited by the planes π 1 and π 2 , which are orthogonal to light vectors s 1 and s 2 . (b) An example where the system (8) admits two solutions in n (marked in red).

Differential Formulations of the PS2 Problem

Taking m = 2 in [START_REF] Sato | Reflectance analysis under solar illumination[END_REF], and adding a Dirichlet boundary condition, a first differential formulation of the PS2 problem is written as two nonlinear PDEs:

               ρ(x, y) -s 1 • ∇u(x, y) + s 1 3 1 + ∇u(x, y) 2 = I 1 (x, y) a.e. (x, y) ∈ Ω ρ(x, y) -s 2 • ∇u(x, y) + s 2 3 1 + ∇u(x, y) 2 = I 2 (x, y) a.e. (x, y) ∈ Ω u(x, y) = g(x, y) ∀(x, y) ∈ ∂Ω (10)
Besides its nonlinearity, the main drawback of Problem [START_REF] Onn | Integrability Disambiguates Surface Recovery in Two-Image Photometric Stereo[END_REF] concerns the need for a boundary condition: the function g(x, y), taken in the space of Lipschitz functions, represents a piece of information which is rarely available. On the other hand, since a common factor ρ(x, y)/ 1 + ∇u(x, y) 2 occurs in both PDEs of [START_REF] Onn | Integrability Disambiguates Surface Recovery in Two-Image Photometric Stereo[END_REF], we can combine them in order to simultaneously eliminate the albedo ρ(x, y) and the nonlinearity (we only have to suppose that ρ(x, y) = 0):

I 2 (x, y) s 1 -I 1 (x, y) s 2 • ∇u(x, y) = I 2 (x, y) s 1 3 -I 1 (x, y) s 2 3 ( 11
)
Considering the same boundary condition as in [START_REF] Onn | Integrability Disambiguates Surface Recovery in Two-Image Photometric Stereo[END_REF], we deduce:

b(x, y) • ∇u(x, y) = f (x, y) a.e. (x, y) ∈ Ω u(x, y) = g(x, y) ∀(x, y) ∈ ∂Ω (12)
where:

b(x, y) = I 2 (x, y) s 1 -I 1 (x, y) s 2 f (x, y) = I 2 (x, y) s 1 3 -I 1 (x, y) s 2 3 ( 13
)
This second differential formulation of the PS2 problem allows us to propagate boundary information g(x, y) across Ω through vector field b(x, y). Indeed, if b(x, y) and f (x, y) are two bounded (but not necessarily continuous) functions defined by ( 13), and if g(x, y) is a Lipschitz function, then ( 12) admits a unique Lipschitz solution u(x, y) [START_REF] Mecca | Uniqueness and Approximation of a Photometric Shape-from-Shading Model[END_REF]. It is worth emphasizing that this implies that the differential formulations of PS2 allow one to reconstruct surfaces that are differentiable almost everywhere. In practice, this means that surfaces having sharp structures could be retrieved. In contrast, using the non-differential approach, the surface must be assumed to be C 1 , so that the normal is defined everywhere. If edges or depth discontinuities are in fact present, they must be adequately handled during the integration stage [START_REF] Durou | Integration of a Normal Field in the Presence of Discontinuities[END_REF].

In order to ensure robustness to noise, propagation schemes used in [START_REF] Mecca | Unambiguous Photometric Stereo Using Two Images[END_REF] may be advantageously replaced by variational methods, recasting the linear PDE [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF] as the following optimization problem:

   min u: Ω→R Ω b(x, y) • ∇u(x, y) -f (x, y) 2 dx dy s.t. u(x, y) = g(x, y), ∀(x, y) ∈ ∂Ω (14)
When m ≥ 3 images are available, it was shown in [START_REF] Quéau | Unbiased photometric stereo for colored surfaces: A variational approach[END_REF] that variational models such as ( 14) may be considered in several more difficult situations such as color PS, PS with pointwise sources or perspective PS, even without a boundary condition. Yet, when m = 2, explicit knowledge of the function g is required in order to ensure that the characteristics are not reconstructed independently. An alternative consists in using anisotropic regularization, in order to "couple" the 3D-reconstructions of the different characteristics, by ensuring smoothness along directions that are not tangent to the characteristics. A suitable variational model, which was introduced by Hernández et al. in [START_REF] Hernández | Overcoming shadows in 3-source photometric stereo[END_REF], is written as follows:

min u: Ω→R Ω b(x, y) • ∇u(x, y) -f (x, y) 2 + α 1 b ⊥ (x, y) • ∇u(x, y) 2 + α 2 b ⊥ (x, y) H(u)(x, y)b ⊥ (x, y) 2 dxdy ( 15
)
where the field b ⊥ is perpendicular to the characteristic curves and H(u) is the Hessian matrix of u. The parameters α 1 and α 2 must be tuned appropriately, in order to ensure that the regularization is sufficient, yet avoiding oversmoothing the solution (cf. Section 5.4).

Example of PS2 Problem

In order to question the consistency between the different formulations of the PS2 problem, let us take the example of a plane surface u(x, y) = x illuminated by light vectors s 1 = [0, 0, 1] and

s 2 = 1 2 [1, 1, √ 2] . If ρ ≡ 1
, the images of this surface using Model (6) are uniform, with graylevels

I 1 (x, y) = 1 √ 2 and I 2 (x, y) = √ 2-1 2 √ 2 .
Let us first solve this PS2 example using the non-differential formulation [START_REF] Jung | One-day outdoor photometric stereo via skylight estimation[END_REF], which is written as:

             n 3 (x, y) = 1 √ 2 n 1 (x, y) + n 2 (x, y) + √ 2 n 3 (x, y) 2 = √ 2 -1 2 √ 2 n 1 (x, y) 2 + n 2 (x, y) 2 + n 3 (x, y) 2 = 1 (16)
This system admits two solutions n

+ = 1 √ 2 [-1, 0, 1] and n -= 1 √ 2 [0, -1, 1] ,
independently of (x, y). If the surface to be reconstructed is supposed to be differentiable everywhere, there are only two acceptable normal fields n(x, y) = n + or n(x, y) = n -. It follows from (4) that there are two possible values for ∇u(x, y):

-1 n + 3 (x, y) n + 1 (x, y) n + 2 (x, y) = 1 0 ; -1 n - 3 (x, y) n - 1 (x, y) n - 2 (x, y) = 0 1 (17)
Imposing u(0, 0) = 0 to fix the integration constants, we finally obtain two solutions u + (x, y) = x and u -(x, y) = y. Fortunately, one of these solutions is the genuine surface. Now, let us write the first differential formulation [START_REF] Onn | Integrability Disambiguates Surface Recovery in Two-Image Photometric Stereo[END_REF], supposing ρ ≡ 1:

         1 1 + ∇u(x, y) 2 = 1 √ 2 -∂ x u(x, y) -∂ y u(x, y) + √ 2 2 1 + ∇u(x, y) 2 = √ 2 -1 2 √ 2 (18)
Note that no boundary condition is available. System (18) is equivalent to:

∂ x u(x, y) + ∂ y u(x, y) = 1 ∇u(x, y) 2 = 1 (19)
By replacing ∇u(x, y) 2 with ∂ x u(x, y) 2 +∂ y u(x, y) 2 , we quickly find that (19) has the same solutions [START_REF] Yang | Two image photometric stereo method[END_REF] in ∇u(x, y) as ( 16). Both formulations are thus consistent, in the sense that they provide the same global solutions. It follows from the definitions (13) that:

b(x, y) = - 1 2 √ 2 [1, 1] ; f (x, y) = - 1 2 √ 2 (20)
In the absence of a boundary condition, the second differential formulation [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF] reduces to:

∂ x u(x, y) + ∂ y u(x, y) = 1 (21)
Unsurprisingly, this PDE is the first equation of ( 19), which admits an infinity of solutions. For instance, all functions u(x, y) = a x + (1 -a) y + w(x -y), for any a ∈ R and any differentiable function w : R → R, are solutions to (21). How can it be explained that this formulation does not lead to the same conclusion as the others? In fact, without a boundary condition, differential formulation ( 12) is a necessary but insufficient condition. For better constraint, a boundary condition is required (but rarely available).

To conclude, there is obviously no inconsistency between the non-differential formulation and the first differential formulation of the PS2 problem. Hence, both these formulations should allow us to predict the same number of global solutions. Let us examine in more detail this issue.

PS2 Problem: Predicting the Number of Solutions

At any twice-lit point (x, y) ∈ Ω, we know from Section 3.1 that the PS2 problem admits one or two solutions in n(x, y). Of course, in order to predict the number of normal fields, the points where the normal can be determined unambiguously are of primary importance. Such singular points have been studied in detail to solve the SFS problem [START_REF] Ikeuchi | Numerical Shape from Shading and Occluding Boundaries[END_REF].

Singular Points

As already observed in Section 3.1, the first situation where normal uniqueness can be proved is reached when both solutions n + (x, y) and n -(x, y) to the problem (8) coincide. In this case, the unique solution is inside plane π. The set S R of such normals is a geodesic on G. It is the intersection between S and π (see Fig. 2-a). A second type of singular point occurs when one solution to (8) points towards the viewer (S G set in Fig. 2-b), while the other one points away (S Y set in Fig. 2-b). In this case, the ambiguities are easy eliminated by choosing the "visible" normal. Let Ω R and Ω G be the singular points sets where the normal belongs to S R or S G , respectively. Since a PS2 problem comes down to a pair of SFS problems (with the constraint that the camera pose is unique), one could wonder why the singular points of each SFS problem are not taken into account. These points are such that n = s 1 in the first image, or n = s 2 in the second. However, both these values of n are inside S R , since they support plane π. This shows us that the singular points of the PS2 problem include those of each SFS subproblem.

s 1 s 2 S R E z x y S G S Y s 1 s 2 z y x E (a) (b)
The normal field is continuous if the surface is supposed to be C 1 . Under such an assumption, is it possible to propagate the knowledge of the normal in a singular point to its non-singular neighbours? The answer to this question depends on which type of singular points is referred to. Let P R ∈ Ω R be a singular point of the first type i.e., one whose normal is inside S R . The normal in a non-singular neighbour PR of P R can lie on both sides of S R . The subsets of S which are above and below the geodesic S R are respectively called S U and S B , see Fig. 3. In other words, there is a remaining ambiguity on the normal in PR . Now, let P G ∈ Ω G be a singular point of the second type, whose normal is inside S G . In any non-singular neighbour PG of P G , we can infer from the normal field continuity that the normal is inside S U . To conclude, the normal is unambiguously known in all points connected to Ω G inside Ω\Ω R , which can be considered as supplementary singular points. The sets of singular points Ω R and Ω G , which are numerically estimated from the estimation of n + and n -described in Section 5, are superimposed to these images. In each of the five connected parts of Ω\Ω R , there are two solutions: one is contained in S U , the other in S B . This gives rise to 2 5 = 32 continuous normal fields. Nevertheless, a difference between our count and that put forward by Onn and Bruckstein [START_REF] Onn | Integrability Disambiguates Surface Recovery in Two-Image Photometric Stereo[END_REF] comes from the detection of Ω G : since the region Ω\Ω R containing Ω G is determined, only the other four regions are ambiguous, leading eventually to 2 4 = 16 possible normal fields.

s 1 s 2 S B S R E z x y s 2 s 1 S G S U S R S Y y x E (a) (b)
The solution u can be calculated over each connected singular points set, up to a constant of integration. Knowing that differential formulation [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF] needs a boundary condition to be well-posed, this is a simple way to construct one. Since ∂Ω is connected to Ω G inside Ω\Ω R (see Figs. 4-b and4-c), the solution can be univocally calculated along ∂Ω. According to [START_REF] Mecca | Uniqueness and Approximation of a Photometric Shape-from-Shading Model[END_REF], this allows us to predict a unique solution, which is of course the original surface shown in Fig. 4-a. Based on this example, we could conclude that the non-differential and differential approaches to the PS2 resolution are complementary, in order to predict the number of solutions. Namely, the boundary condition required by the latter is provided by the former.

A Possible Remaining Ambiguity

The conclusion of the previous section does not always hold true. For instance, in a case such as the example of Section 3.3, there is no singular point (all the points have the same normal). Let us show another counterexample originally exhibited by Kozera in [START_REF] Kozera | Existence and Uniqueness in Photometric Stereo[END_REF]. The surface represented by equation z = x y, with uniform albedo ρ ≡ 1, illuminated by light vectors s 1 = [s, s, c] and s 2 = [-s, -s, c] , where s = √ 2/2 sin φ and c = cos φ, for a given φ ∈ ]0, π/2[, is characterized by the following graylevel functions:

         I 1 (x, y) = s (-x -y) + c 1 + x 2 + y 2 I 2 (x, y) = -s (-x -y) + c 1 + x 2 + y 2 (22)
The two first equations of the non-differential problem (8) are rewritten as:

         s [n 1 (x, y) + n 2 (x, y)] + c n 3 (x, y) = s (-x -y) + c 1 + x 2 + y 2 -s [n 1 (x, y) + n 2 (x, y)] + c n 3 (x, y) = -s (-x -y) + c 1 + x 2 + y 2 (23)
Since by definition, s and c are nonzero, System (23) is equivalent to:

       n 1 (x, y) + n 2 (x, y) = -x -y 1 + x 2 + y 2 n 3 (x, y) = 1 1 + x 2 + y 2 (24)
Using (24), the third equation of ( 8) can be rewritten as:

n 1 (x, y) 2 + x + y 1 + x 2 + y 2 n 1 (x, y) + x y 1 + x 2 + y 2 = 0 ( 25
)
It is easy to check that this second-order equation always admits two real solutions in n 1 (x, y), which come down to a unique solution when y = x.

These solutions give rise to two possible normals at each point (x, y) ∈ Ω:

n + (x, y) = 1 1 + x 2 + y 2   -y -x 1   ; n -(x, y) = 1 1 + x 2 + y 2   -x -y 1   (26)
We deduce from ( 26) and ( 4) two possible values for ∇u(x, y):

-1 n + 3 (x, y) n + 1 (x, y) n + 2 (x, y) = y x ; -1 n - 3 (x, y) n - 1 (x, y) n - 2 (x, y) = x y (27)
Both these vector fields are easily integrated, which provides us with two solutions u + (x, y) = x y and u -(x, y) = (x 2 + y 2 )/2, up to two additive constants. That is to say, there is a remaining ambiguity.

On the other hand, Ω G is empty, but it is easily deduced from ( 26) that Ω R is the straight line y = x. The solution can therefore be calculated along this line, up to a constant. Using a similar rationale as in Section 4.2, we should conclude that the solution is unique, which would contradict the previous result. This contradiction is easily explained: the prediction of [START_REF] Mecca | Uniqueness and Approximation of a Photometric Shape-from-Shading Model[END_REF] holds true only if the solution is known on a curve which is not a characteristic. As we will see in the next section, this condition is precisely not valid in this case.

Using the Characteristics to Predict the Number of Solutions

In the previous example, each vector b(x, y) defined in Eq. ( 13) is parallel to [START_REF]Shape from Shading[END_REF][START_REF]Shape from Shading[END_REF] . The characteristics are thus the straight lines represented by equations y = x + η, η ∈ R, including the set Ω R . Depth u can be univocally calculated along each characteristic, up to a constant. This uniqueness result is not contradicted by the previous two-fold ambiguity since, for any η ∈ R:

u + (x, x + η) -u -(x, x + η) = x(x + η) - x 2 + (x + η) 2 2 = - η 2 2 ( 28
)
is independent from x.

If u(x, y) is known at one point of each characteristic, we have a better understanding why Problem (12) has a unique solution. Following this rationale, all functions of the following form seem to be solutions to the previous example:

u(x, y) = u + (x, y) + v(y -x) (29)
provided that v is a scalar function such that v(y -x) is constant along each characteristic. In fact, any function v is not acceptable because, as already noted, differential formulation ( 12) is a necessary but insufficient condition, in the absence of boundary condition. From (29), we deduce:

∇u(x, y) = y -v (y -x) x + v (y -x) (30)
Eq. ( 7) tells us that the surface z = u(x, y) of albedo ρ ≡ 1 is a solution to the previous example only if:

         -s [y -v (y -x)] -s [x + v (y -x)] + c 1 + [y -v (y -x)] 2 + [x + v (y -x)] 2 = I 1 (x, y) s [y -v (y -x)] + s [x + v (y -x)] + c 1 + [y -v (y -x)] 2 + [x + v (y -x)] 2 = I 2 (x, y) (31)
Using [START_REF] Ackermann | Photometric Stereo for Outdoor Webcams[END_REF], we easily find the only two solutions in v (y -x) to System (31):

v 1 (y -x) = 0 v 2 (y -x) = y -x =⇒    v 1 (y -x) = K 1 , K 1 ∈ R v 2 (y -x) = (y -x) 2 2 + K 2 , K 2 ∈ R (32)
Plugging ( 32) into (29), we eventually obtain the two following solutions:

   u 1 (x, y) = u + (x, y) + K 1 u 2 (x, y) = u + (x, y) + (y -x) 2 2 + K 2 = u -(x, y) + K 2 (33)
This result confirms that there are only two analytical solutions u + and u -, up to the constants K 1 and K 2 .

Integrability Constraint

Even if vector field [p, q] = [-n 1 /n 3 , -n 2 /n 3 ] is easily calculated from a normal field, there is no guarantee that this vector field is integrable i.e., that it satisfies the integrability constraint [START_REF] Horn | The Variational Approach to Shape From Shading[END_REF]:

∂p ∂y = ∂q ∂x ( 34
)
whereas this is required, if the surface is supposed to be at least C 2 , in order to be sure that equation ∇u = [p, q] has a solution in u. Onn and Bruckstein note in [START_REF] Onn | Integrability Disambiguates Surface Recovery in Two-Image Photometric Stereo[END_REF] that "most of the time", in each connected part P of Ω\Ω R , one of the two possible normal fields may be discarded since it is not integrable.

To decide, they rely on a criterion deduced from (34):

(x,y)∈P ∂p ∂y (x, y) - ∂q ∂x (x, y) 2 dx dy = 0 ( 35
)
They also characterize the "rare cases" where ( 35) is satisfied by more than one normal field, in which case the PS2 problem admits several solutions. In fact, the examples of Sections 3.3 and 4.3 are such "rare cases".

Let us now return to the example of Fig. 4. Using the non-differential formulation, we found 2 4 = 16 possible normal fields. A more complex rationale based on the topology of sets Ω R and Ω G allowed us to predict a unique solution. Consequently, among the sixteen normal fields, only one is integrable. However, a prediction based on topology is difficult to extend to the discrete framework, since the notion of continuity will be lost.

Indeed, let us calculate the sets Ω R and Ω G when the images are quantized. It is clear from Fig. 5 that the topology of Ω R is very sensitive to quantization noise. It is thus no longer possible to predict the number of solutions using the rationale of Section 4.2 from such a fragmented set Ω R . We show in the next section that it is possible to efficiently find the most integrable normal field, and hence to eliminate the ambiguities of the PS2 problem, without knowledge of a boundary condition [START_REF] Mecca | Unambiguous Photometric Stereo Using Two Images[END_REF] nor parameter tuning [START_REF] Hernández | Overcoming shadows in 3-source photometric stereo[END_REF].

Photometric Stereo Using Two Images: A Numerical Resolution

We know from Section 3.1 that the PS2 problem admits at most two solutions in n at each point (x, y) ∈ Ω (see Fig. 6). We can thus a priori construct at most 2 |Ω| different discrete normal fields, where |Ω| denotes the number of pixels inside Ω. Yet, except in some rare cases [START_REF] Onn | Integrability Disambiguates Surface Recovery in Two-Image Photometric Stereo[END_REF], only one of these candidates is integrable. Let us show how to efficiently find this most integrable normal field. 42)), using the quantized dataset presented in Fig. 5. Bottom: 3D-shapes obtained integrating the normals [START_REF] Durou | Integration of a Normal Field in the Presence of Discontinuities[END_REF]. The normal fields n + and n -are only two solutions among 2 |Ω| , since any combination of both these normal fields is also plausible.

n n + n - u u + u -

Estimating the Candidate Normal Fields

Problem (8) has two solutions n + and n -at each point. Let us first show how to express these solutions using a purely algebraic method.

At each point (x, y) ∈ Ω, the two first equations of Problem ( 8) write (the dependencies in (x, y) are omitted):

s 1 s 2 S ∈R 2×3 n = I 1 ρ I 2 ρ i∈R 2 (36)
With a view to solving Eq. ( 36) in the least-squares sense in order to handle quantization noise, we reformulate problem (8) as follows:

min n S n -i 2 s.t. n = 1 (37)
The singular value decomposition (SVD) of S can be written as:

S = UΣV = σ 1 u 1 v 1 + σ 2 u 2 v 2 (38)
In [START_REF] Vogel | Making Shape from Shading Work for Real-World Images[END_REF], σ 1 ≥ σ 2 > 0 are the pair of strictly positive singular values of S (S has rank 2, since s 1 and s 2 are non-collinear), u 1 and u 2 are orthonormal vectors of R 2 , and v 1 and v 2 are orthonormal vectors of R 3 . It can be shown (Theorem 5.5.1 in [START_REF] Golub | Matrix computations[END_REF]) that the minimum-norm solution of the least-squares problem resulting from (36):

   n 0 = argmin n∈E n s.t. E = n ∈ R 3 , S n -i 2 ≤ S m -i 2 ∀m ∈ R 3 (39)
is written as:

n 0 = u 1 i σ 1 v 1 + u 2 i σ 2 v 2 (40)
Let us now solve Problem [START_REF] Courteille | Towards shape from shading under realistic photographic conditions[END_REF]. Three cases can occur:

1. If n 0 = 1, then n + = n -= n 0 is the only solution to [START_REF] Courteille | Towards shape from shading under realistic photographic conditions[END_REF], and the objective is null. This corresponds to the singular points in Ω R , cf. Section 4.1. 2. If n 0 < 1, then neglecting the unit-length constraint, there exists a set of solutions with null objective to the optimization problem [START_REF] Courteille | Towards shape from shading under realistic photographic conditions[END_REF], which is written as:

n 0 + Rv 3 ( 41
)
where v 3 is a unit-length vector of the kernel of S (this kernel has dimension 1, according to the rank theorem). Two of these solutions have unit-length, which are the only two solutions to Problem (37)2 :

n + = n 0 + 1 -n 0 2 v 3 n -= n 0 -1 -n 0 2 v 3 (42)
3. Finally, if n 0 > 1, there is no solution with null objective to problem [START_REF] Courteille | Towards shape from shading under realistic photographic conditions[END_REF]. It is tempting to choose the approximate solution n + = n -= n 0 n 0 , but this is not the real solution to [START_REF] Courteille | Towards shape from shading under realistic photographic conditions[END_REF]. It is proven in [START_REF] Golub | Matrix computations[END_REF] that this solution writes as:

n(λ 0 ) = σ 1 u 1 i σ 2 1 + λ 0 v 1 + σ 2 u 2 i σ 2 2 + λ 0 v 2 (43)
which reduces to n 0 when λ 0 = 0, yet λ 0 is in fact the unique positive solution in λ to the following secular equation:

σ 1 u 1 i σ 2 1 + λ 2 + σ 2 u 2 i σ 2 2 + λ 2 -1 = 0 (44)
Eq. ( 44) could be rewritten under the form of an algebraic equation of degree 4 in λ, and solved using the Ferrari-Cardan formulae. Instead, we used in our implementation a Newton method. However, it should be reminded that an equation of this type has to be solved at each pixel such that n 0 > 1. If the targetted application has real-time requirements [START_REF] Hernández | Nonrigid Photometric Stereo with Colored Lights[END_REF], n 0 / n 0 may be preferred as a fast approximation of the solution to [START_REF] Courteille | Towards shape from shading under realistic photographic conditions[END_REF], although it is not the exact solution.

Disambiguating the Problem by Graph Cut

The case n 0 < 1 being by far the most frequent, we can actually build almost 2 |Ω| normal fields, which are all solutions to Problem [START_REF] Jung | One-day outdoor photometric stereo via skylight estimation[END_REF]. In [START_REF] Woo | 3D shape restoration using sparse representation and separation of illumination effects[END_REF], it is proposed to better constrain the problem assuming that the normals are distributed according to a Laplace law, but this assumption is hard to justify. As discussed in Section 4.5, we would rather advise using the integrability constraint of the normal field, which is much less restrictive (the surface is simply assumed C 2 , at least piecewise). We will now describe the practical resolution of the problem (8) based on this constraint.

Finding the "most integrable" normal field amounts to finding the one such that the integrability constraint [START_REF] Golub | Matrix computations[END_REF] is "best" approximated over Ω. This can be formulated as the variational problem [START_REF] Woo | 3D shape restoration using sparse representation and separation of illumination effects[END_REF]. Yet, such a variational problem does not account for the fact that we know explicitly the 2 |Ω| possible normal fields. An exhaustive search of the most integrable normal field could be preferred, but it may be computationally infeasible on large data. Instead, we suggest to recast the task of finding the most integrable normal field as a binary labelling problem whose solution can be found efficiently, by means of the graph cut algorithm [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF].

Let us attribute to each pixel a label l ∈ {+, -} indicating the normal n + or n -, and let us denote by [p l , q l ] = [-n l 1 /n l 3 , -n l 2 /n l 3 ] the corresponding discrete approximation of the surface gradient ∇u. The "optimal" labelling l : Ω → {+, -} is the one which makes the normal field the "most integrable". Thus, we now have to solve the following discrete version of the variational problem [START_REF] Woo | 3D shape restoration using sparse representation and separation of illumination effects[END_REF] over Ω:

min l (x,y)∈Ω ∂p l ∂y (x, y) - ∂q l ∂x (x, y) 2 (45)
In order to discretize the space derivatives of Problem (45) using finite differences of order 1, we should consider the four possible cliques families of order 3:

         C 3 1 = {(x, y), (x -1, y), (x, y -1)} ∈ Ω 3 C 3 2 = {(x, y), (x + 1, y), (x, y -1)} ∈ Ω 3 C 3 3 = {(x, y), (x -1, y), (x, y + 1)} ∈ Ω 3 C 3 4 = {(x, y), (x + 1, y), (x, y + 1)} ∈ Ω 3 (46)
If we denote (l 1 , l 2 , l 3 ) the labels of each clique's three pixels, in the same order they are defined in [START_REF] Papadhimitri | A New Perspective on Uncalibrated Photometric Stereo[END_REF], Problem ( 45) can be rewritten as:

min l c 3 ∈C 3 V int c 3 (l 1 , l 2 , l 3 ) ( 47
)
where C 3 = ∪ 4 i=1 C 3 i , and potential V int c 3 (l 1 , l 2 , l 3 ) gives the local integrability for the current clique c 3 and the current labeling l. For example, if c 3 1 ∈ C 3 1 :

V int c 3 1 (l 1 , l 2 , l 3 ) = p l 1 (x, y) -p l 3 (x, y -1) -q l 1 (x, y) -q l 2 (x -1, y) 2 (48)
Problem ( 47) is a labeling problem where the local potential depends on the current pixel and on two of its neighbours. Such combinatorial optimization problems have been studied in [START_REF] Kolmogorov | What Energy Functions Can Be Minimized via Graph Cuts?[END_REF], where it has been proven that the graph cut algorithm [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF] can be used to minimize the energy c 3 ∈C 3 V int c 3 (l 1 , l 2 , l 3 ), provided that its regularity (sub-modularity) is ensured, which means here:

     V int c 3 (+, +, l 3 ) + V int c 3 (-, -, l 3 ) ≤ V int c 3 (+, -, l 3 ) + V int c 3 (-, +, l 3 ) V int c 3 (+, l 2 , +) + V int c 3 (-, l 2 , -) ≤ V int c 3 (+, l 2 , -) + V int c 3 (-, l 2 , +) V int c 3 (l 1 , +, +) + V int c 3 (l 1 , -, -) ≤ V int c 3 (l 1 , +, -) + V int c 3 (l 1 , -, +) (49)
for any c 3 ∈ C 3 and any (l 1 , l 2 , l 3 ) ∈ {+, -} 3 . Of course, these inequalities have no reason to be satisfied.

Ensuring the Regularity Condition

To ensure the regularity condition, we modify the problem (47) by introducing a regularization term of the Ising type:

min l c 3 ∈C 3 V int c 3 (l 1 , l 2 , l 3 ) + c 2 ∈C 2 V Ising c 2 (l 1 , l 2 ) ( 50
)
where C 2 = ∪ 4 i=1 C 2 i gathers the four following sets of cliques of order 2:

         C 2 1 = {(x, y), (x -1, y)} ∈ Ω 2 C 2 2 = {(x, y), (x, y -1)} ∈ Ω 2 C 2 3 = {(x, y), (x -1, y -1)} ∈ Ω 2 C 2 4 = {(x, y), (x + 1, y -1)} ∈ Ω 2 (51) and V Ising c 2
(l 1 , l 2 ) is defined as follows, where (l 1 , l 2 ) denote the labels of the two pixels of each clique c 2 ∈ C 2 , in the same order they are defined in (51):

V Ising c 2 (l 1 , l 2 ) = β c 2 δ(l 1 = l 2 ) ( 52
)
In Eq. ( 52), β c 2 is a positive or zero local coefficient i.e., such a coefficient must be fixed for each c 2 ∈ C 2 . Enforcing the regularity condition for (50), we deduce a lower bound for each β c 2 . For example, if c 2 3 = {(x, y), (x -1, y -1)} ∈ C 2 3 :

β c 2 3 ≥ 1 2 max 0, max (l 3 ,l 4 )∈{+,-} 2 ∆V int c 3 3 (l 3 ), ∆V int c 3 2 (l 4 ) ( 53
)
where c3 3 = c 2 3 ∪ (x, y -1), c 3 2 = c 2 3 ∪ (x -1, y), and (l 3 , l 4 ) denote the labels of (x, y -1) and (x -1, y), respectively. Finally, for j ∈ {2, 3} and k ∈ {+, -}:

∆V int c 3 j (k) = V int c 3 j (k, +, +) + V int c 3 j (k, -, -) -V int c 3 j (k, +, -) -V int c 3 j (k, -, +) (54)
From a Markovian point of view, our approach consists in using nonstationary Ising models. The use of such models corresponds to a piecewise uniform prior on the labeling l. We know from Section 4 that there exist a finite number of connected areas over which the solutions n + and n -are different (these areas are bounded by the set Ω R ). If these normal fields are estimated as indicated in Section 5.1, they will be continuous inside each area, so the optimal labeling should change only along their boundaries. The Ising prior is thus consistent with the theoretical analysis conducted in Section 4.

Yet, in order not to bias the 3D-reconstruction which should rather be guided by integrability, this prior should have as little influence as possible: according to inequalities like (53), it is possible to predict, for each clique c 2 ∈ C 2 , the smallest value of coefficient β c 2 to ensure the regularity condition. These coefficients should therefore not be considered as parameters, which is an advantage. This allows us to limit the energy regularization, only in order to ensure regularity in the sense of Kolmogorov [START_REF] Kolmogorov | What Energy Functions Can Be Minimized via Graph Cuts?[END_REF], thus avoiding oversmoothing. When coefficients β c 2 are not fixed to their minimum values, the labeling is oversmoothed, which biases the 3D-reconstruction (see Fig. 7). Our approach thus avoids the difficulty of tuning the regularization parameters, which is a problem with existing methods (cf. Figs. 8 and9).

An Efficient Method for Solving the PS2 Problem

To sum up, the method of resolution of the PS2 problem that we recommend comprises two stages:

1. The calculation of normal fields n + and n -as indicated in Section 5.1.

The disambiguation of the problem using the integrability criterion, an

8-connected Ising model, the minimum values of the coefficients β c 2 like for instance (53), and the graph cut algorithm.

Both steps can be efficiently conducted: the initial solutions n + and n -are explicit, and optimization by graph cut is very efficient 3 . The 3D-reconstruction shown in the first row of Fig. 7 is obtained applying this method to the example of Fig. 6. The normal field is very similar to that of the ground truth: when the albedo is known and in the absence of shadow, our method is able to recover almost exactly the genuine normals.

In addition to being parameter-free and not requiring a boundary condition, Figs. 8 and9 show us that the proposed method is able to confine estimation errors, which may be useful if a shadow is present in one (or both) of the images. This is in contrast with existing methods, for which an estimation error in one pixel may propagate to its neighbors. 5 was considered while increasing the zenithal angles φ for the lightings, in order to create stronger shadowing effects. Our method confines the errors in the shadow areas, while requiring neither a boundary condition nor parameter tuning (α 1 and α 2 are defined in Eq. ( 15)).

[31] + Dirichlet [START_REF] Hernández | Overcoming shadows in 3-source photometric stereo[END_REF] (α 1 = 0.15, α 2 = 1)

[9] (α 1 = 0.01, α 2 = 0.3) Proposed Figure 9: 3D-reconstructions corresponding to the second row in Fig. 8.

Yet, our approach has one major weakness with respect to existing methods. Indeed, we assume that the albedo is known. In contrast, since existing methods [START_REF] Hernández | Overcoming shadows in 3-source photometric stereo[END_REF][START_REF] Quéau | Unbiased photometric stereo for colored surfaces: A variational approach[END_REF] are based on image ratios, knowledge of the albedo is not required. As a consequence, their methods are robust to unpredicted albedo variations, while ours is not. As it is illustrated in Figs. 10 and11, this represents an important limitation of our method. To conclude, our method yields overall more accurate 3D-reconstructions than existing ones, and is more flexible (neither need for boundary condition nor parameter tuning), but only when the albedo is known. This restriction is also recurrent in shapefrom-shading, where albedo estimation must be carried out beforehand using interpolation techniques [START_REF] Courteille | Towards shape from shading under realistic photographic conditions[END_REF][START_REF] Vogel | Making Shape from Shading Work for Real-World Images[END_REF], or within the 3D-reconstruction process by introducing priors on the albedo and the shape [START_REF] Barron | High-Frequency Shape and Albedo from Shading using Natural Image Statistics[END_REF].

For sake of completeness, well-posedness of the PS2 problem with unknown albedo is briefly discussed hereafter. The 256 × 256 images were created with ρ everywhere equal to 1 except on a diagonal band with width d ρ pixels where it was set to 0.9. Assuming (wrongly) that ρ ≡ 1 induces a bias with our method, while existing methods are independent from the albedo.

[31] + Dirichlet [START_REF] Hernández | Overcoming shadows in 3-source photometric stereo[END_REF] (α 1 = 0.15, α 2 = 1)

[9] (α 1 = 0.01, α 2 = 0.3) Proposed

Figure 11: 3D-reconstructions corresponding to the second row in Fig. 10.

PS2 Problem with Unknown Albedo

At this stage, it is interesting to quote a variant of the PS2 problem, where the surface to be reconstructed has an unknown albedo. The presence of a supplementary unknown ρ(x, y) at each point (x, y) ∈ Ω could drastically complicate the problem. A main feature of the PS3 problem is that the albedo can be univocally estimated [6]. This is not the case for SFS since, even if the albedo is known, the problem is usually ill-posed [START_REF] Oliensis | Uniqueness in Shape from Shading[END_REF]. We will see that the PS2 problem with unknown albedo is more similar to SFS than to PS3.

Let us first give a geometric interpretation of [START_REF] Mecca | Unambiguous Photometric Stereo Using Two Images[END_REF]. Since s = s , s 3 and n(x, y) is parallel to -∇u(x, y) , 1 , this equation is rewritten as:

I 2 (x, y) s 1 -I 1 (x, y) s 2 • n(x, y) = 0 (55)
Eq. ( 55) can be directly derived from the two first equations of [START_REF] Jung | One-day outdoor photometric stereo via skylight estimation[END_REF]. It means that n(x, y) lies within a plane χ(x, y) which is orthogonal to the vector I 2 (x, y) s 1 -I 1 (x, y) s 2 . Naturally, χ(x, y) is orthogonal to the plane π as well, which is supported by the vectors s 1 and s 2 . It is noteworthy that Eq. ( 55) holds true for any albedo value ρ(x, y), because Eq. ( 11) has been derived from Problem [START_REF] Onn | Integrability Disambiguates Surface Recovery in Two-Image Photometric Stereo[END_REF] by elimination of ρ(x, y). We know from Section 3.1 that, for any known albedo value ρ(x, y), Problem (8) has two solutions n + (x, y) and n -(x, y) which are symmetric with respect to π. The geometric interpretation of Eq. ( 55) means that for a fixed pair of graylevels (I 1 (x, y), I 2 (x, y)), different values of ρ(x, y) will provide different values of n + (x, y) and n -(x, y) which all lie within χ(x, y). Fig. 12 shows us that n + (x, y) and n -(x, y) move towards plane π as ρ(x, y) decreases from 1 towards a limit value ρ inf (x, y) which corresponds to the limit case n + (x, y) = n -(x, y). Hence, at each (x, y), a range of values [ρ inf (x, y), 1] are feasible for ρ(x, y). If s 1 = s 2 , it can be shown that this limit is:

ρ inf (x, y) = I 1 (x, y) 2 + I 2 (x, y) 2 -2 I 1 (x, y) I 2 (x, y) (s 1 • s 2 ) 1 -(s 1 • s 2 ) 2 (56)
Hence, the PS2 problem with unknown albedo is ill-posed. However, it is worth underlining that, when a boundary condition is available, Problem [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF] can be solved without any knowledge of ρ(x, y). Moreover, ρ(x, y) can be a posteriori calculated in this case, using any of the PDEs of [START_REF] Onn | Integrability Disambiguates Surface Recovery in Two-Image Photometric Stereo[END_REF]. (x,y) move towards π as ρ(x, y) decreases from 1 towards a minimum value ρ inf (x, y) which corresponds to the limit case n + (x, y) = n -(x, y).

As an example, let us reconsider the example of Section 3.3, under the assumption that ρ(x, y) is unknown. The non-differential problem (16) becomes:

             ρ(x, y) n 3 (x, y) = 1 √ 2 ρ(x, y) n 1 (x, y) + n 2 (x, y) + √ 2 n 3 (x, y) 2 = √ 2 -1 2 √ 2 n 1 (x, y) 2 + n 2 (x, y) 2 + n 3 (x, y) 2 = 1 (57)
For any given value of ρ(x, y), the two first equations of (57), which can be rewritten as n 3 (x, y) = 1 √ 2 ρ(x,y) and n 1 (x, y) + n 2 (x, y) = -1 √ 2 ρ(x,y) , admit an infinity of solutions:

n(x, y) = 1 ρ(x, y)    1 √ 2   -1 0 1   + t   -1 1 0      (58)
which depend on a real parameter t. Replacing the expression (58) of n(x, y) in the third equation of (57), we find after some algebra the following result:

n(x, y) = 1 2 √ 2 ρ(x, y)  -1 - 4 ρ(x, y) 2 -3 -1 + 4 ρ(x, y) 2 -3 2   (59)
where = ±1. In (59), the albedo value ρ(x, y) can be arbitrarily chosen, provided that 4 ρ(x, y) 2 -3 ≥ 0. This means that ρ(x, y) ≥ √ 3/2, which is actually the limit value ρ inf (x, y) given in (56). Clearly, this problem is ill-posed since at each point (x, y) ∈ Ω, for each value ρ(x, y) ∈ [ √ 3/2, 1], there are two possible normals, due to the two possible values of in (59).

To conclude, with the non-differential approach that we follow in this article, the PS2 problem can be solved unambiguously only if the albedo is known beforehand. Nevertheless, let us now show that solving the PS2 problem in an efficient way is not a purely formal challenge.

An Application: Improving Three-source Photometric Stereo

We have seen in Section 5 that the accuracy of PS2 strongly depends on the presence of shadows. The simplest way to ensure robustness to shadows is to consider a third image, i.e. the PS3 problem. By placing the lights appropriately, one can ensure that each surface point is lit in at least two out of the images, and resort to a combination of the PS2 and the PS3 techniques. This can improve a lot the accuracy of real-time PS based on color photometric stereo [START_REF] Hernández | Nonrigid Photometric Stereo with Colored Lights[END_REF].

The Recurrent Problem of Shadows

The three photographs of the first row of Fig. 13, which are available on the web4 , show a plaster bust of Beethoven illuminated by three non-coplanar, parallel and uniform light beams. Since the light vectors are provided, these real data are particularly well adapted to the PS3 technique (see Section 2.2). In addition, let us remark that no boundary condition is available.

Solving at each point (x, y) ∈ Ω a linear system of type (6), then integrating the normals using [START_REF] Durou | Integration of a Normal Field in the Presence of Discontinuities[END_REF], we indeed obtain a "satisfactory" 3D-shape (see Fig. 14), but this is contradicted by the estimated albedo (see Fig. 13-d), which should be uniform since the material is homogeneous. We also note that the points where the albedo estimate is biased lie inside the shadows.

The problem of dealing with shadows in PS is well-known. Because they constitute an unavoidable departure from the Lambertian model (2), shadow graylevels are usually considered as outliers. In this view, most contributions [24, 25,26,[START_REF] Verbiest | Photometric stereo with coherent outlier handling and confidence estimation[END_REF][START_REF] Wu | Robust Photometric Stereo via Low-Rank Matrix Completion and Recovery[END_REF][START_REF] Ikehata | Photometric Stereo Using Sparse Bayesian Regression for General Diffuse Surfaces[END_REF] assume that m > 3 images are available, which is not the case here. Only [START_REF] Hernández | Overcoming shadows in 3-source photometric stereo[END_REF] considers the case m = 3, and thus "two-source photometric stereo [. . . ] in the presence of shadows".

Based on a simple shadow detection (we used the same graph cut-based approach as in [START_REF] Hernández | Shadows in Three-Source Photometric Stereo[END_REF]), Fig. 13-e shows how Ω can be split into four subsets: Ω 3 , Ω 1,2 2 , Ω 1,3 2 and Ω 2,3 2 , with straightforward notations (Ω 1,3 2 is empty in this example). A PS3/PS2 combination can then be considered: the PS2 technique can be used over the pixels lit in only two images, and the PS3 technique elsewhere. The PS3 solution is also used for the points shadowed in more than one image, for stability reasons, yet it would probably be possible to define three other subsets Ω 1 1 , Ω 2 1 and Ω 3 1 and resort to SfS over these sets. To apply the method described in Section 5.4 to the PS3/PS2 combination, we first need to estimate the albedo. Knowing that it is uniform since the material is homogeneous, this can be carried out by evaluating the histogram peak of the estimated albedo inside the set Ω 3 , which can be considered as the "real" albedo: we obtained the value ρ = 0.74.

As shown in Fig. 14, using only two images yields biased results because PS2 is not robust to shadows. Yet, as shown in Fig. 15, combining the PS3 and PS2 techniques can improve a lot the accuracy of the 3D-reconstruction. In addition, it seems that using the proposed PS2 framework inside this combination improves the results, in comparison with state-of-the-art. In order to quantitatively assess the accuracy of the proposed PS2 framework, we extracted three images from the dataset [START_REF] Shi | A Benchmark Dataset and Evaluation for Non-Lambertian and Uncalibrated Photometric Stereo[END_REF], which contains 96 images in total, and for which ground truth normals are available. Our PS2 method was applied only over a shadow area, indicated in red, where the albedo approximately uniform. The proposed method is able to handle shadows even when they are located near the boundaries, which is a known failure case of the differential approach in the absence of a boundary condition [START_REF] Hernández | Overcoming shadows in 3-source photometric stereo[END_REF].

Color Photometric Stereo

Hernández et al. argue the following in [START_REF] Hernández | Overcoming shadows in 3-source photometric stereo[END_REF]: "Using photometric stereo on just three images may seem like an unreasonably hard restriction. There is, however, a particular situation when only three images are available. This technique is known as color photometric stereo". Indeed, the most straightforward application of PS2 is color photometric stereo, a technique where three light sources with different colors and positions are used to simultaneously provide three (graylevel) images of the surface under three different illuminations.

This idea is due to Kontsevich et al., who show in [START_REF] Kontsevich | Reconstruction of shape from shading in color images[END_REF] how to reconstruct the 3D-shape of a white painted scene, illuminated by m = 3 color light sources. Indeed, the number of channels of a standard color image is three. Considering each channel as a graylevel image, a single RGB image is enough to apply the PS3 technique. A deformable scene such as a face can therefore be reconstructed, even if the person is not standing still: Hernández et al. show in [START_REF] Hernández | Nonrigid Photometric Stereo with Colored Lights[END_REF] "how multispectral lighting allows one to essentially capture three images (each with a different light direction) in a single snapshot, thus making per-frame photometric reconstruction possible". An example of such an RGB image extracted from a video sequence 5 is shown in Fig. 17. Let us nevertheless point out that the albedo map should be the same in each channel. In practice, this requires that the scene is made-up. Since applying make-up also ensures that the albedo is uniform, it can be estimated by evaluating the histogram peak. We can apply the PS3 technique to the image of Fig. 17-a but, without any specific treatment, the result is biased around the nose (see Fig. 18), since the shadow renders the red channel unusable in this area (see Fig. 17-b).

As in Section 6.1, the histogram of the albedo allows us to estimate the (uniform) albedo, and therefore to use the PS2 framework within the PS3/PS2 combination. In comparison with PS3, the result is greatly improved, even if artifacts are visible at the junction between the Ω 3 and Ω 2 areas. The proposed method even qualitatively outperforms the state of the art [START_REF] Hernández | Overcoming shadows in 3-source photometric stereo[END_REF].

Conclusion and Future Prospects

In this paper, we first compared in Section 3 the non-differential formulation of the PS2 problem with two differential formulations. This comparison led us to provide in Section 4 a theoretical study of PS2. Building upon this study, we chose the non-differential formulation of PS2 for its numerical resolution, which is detailed in Section 5. With this approach, two normals can explain the graylevels at each point. Our main contribution consists in reformulating the 3D-reconstruction problem as a binary labelling problem, which is solved by means of the graph cut algorithm. The criterion to be optimized is integrability, and a non-stationary Ising term must be added to satisfy the regularity condition.

The straightforward application of this contribution is color photometric stereo, which is explored in Section 6. Our method is validated on RGB images extracted from a video sequence. It improves the 3D-reconstruction fineness, while being quite fast. The real-time requirement expressed in [START_REF] Hernández | Nonrigid Photometric Stereo with Colored Lights[END_REF] is thus not compromised. We moreover assert that this approach is more justified than regularization schemes. Regularization should not be used for disambiguating the 3D-reconstruction problem, but rather for smoothing the residual noise. In this view, future extensions of our work may include an additional regularization designed to ensure more robustness. Yet, designing computationally-efficient algorithms for regularizing normal fields is not an easy task, because of the unit-length constraint [START_REF] Zeisl | Variational Regularization and Fusion of Surface Normal Maps[END_REF]. Recent developments on half-quadratic algorithms designed for manifold-valued images [START_REF] Bergmann | Restoration of manifold-valued images by half-quadratic minimization[END_REF] may constitute an interesting option.

Nevertheless, the assumption that the albedo is known greatly limits the applicability of our study. As future work we plan to extend it to the PS2 problem with unknown albedo. A solution to this more complex problem would be, assuming the albedo known and uniform, that is to say ρ 0 ≡ 1, to get a first estimate u 0 of the 3D-shape using the method described in Section 5.1. It seems then possible to solve the problem, using the constraint u ≈ u 0 as an additional regularization, to deduce a new estimate ρ 1 , and then to loop. A natural perspective of this extension is to cope with images taken by a webcam on a sunny day, as in [START_REF] Jung | One-day outdoor photometric stereo via skylight estimation[END_REF], since the sun rays move within a plane and the linear system (6) is not full-rank.

Finally, throughout this article, we made the hypothesis of an orthographic camera. The proposed method of resolution of the PS2 problem could be used as such for a perspective camera, although the constraint of integrability is more complicated [START_REF] Papadhimitri | A New Perspective on Uncalibrated Photometric Stereo[END_REF], as this only serves to disambiguate a combinatorial optimization problem. Obviously, the perspective should be taken into account in the final step, when the normals are integrated [START_REF] Durou | Integration of a Normal Field in the Presence of Discontinuities[END_REF].

Figure 2 :

 2 Figure 2: (a) Red geodesic S R is the intersection between S and π. Each normal pointing to S R is known without ambiguity. (b) This also holds true for each normal pointing to S G , since the other possible normal points towards the non-visible part S Y of S.

Figure 3 :

 3 Figure 3: (a) Set S B is colored in pink. It is the subset of S bounded by equator E and geodesic S R . (b) Gaussian sphere G is seen from the camera point of view (i.e. z direction).Orange set S U is the subset of S located between S G and S R .

4. 2 .

 2 Using the Singular Points to Construct a Boundary Condition As an illustration, let us calculate a pair of images of the smooth surface depicted in Fig. 4-a, supposing ρ ≡ 1. The light directions are given, respectively, by (θ 1 , φ 1 ) = (60 • , 17 • ) and (θ 2 , φ 2 ) = (135 • , 17 • ), which avoids shadows (see Figs. 4-b and 4-c).

Figure 4 :

 4 Figure 4: (b-c) A pair of 256 × 256 synthetic images (stored in 32 bits) of the smooth surface (a), such that all the points are twice-lit, over which Ω R (in red) and Ω G (in green) are superimposed.

Figure 5 :

 5 Figure5: Two same images as in Fig.4, quantized using 256 levels (8 bits). While this is not the case of Ω G (right), Ω R becomes fragmented (left).

Figure 6 :

 6 Figure 6: Top: RGB-encoded normal fields (n is the ground truth, n + and n -are the normal fields estimated by (42)), using the quantized dataset presented in Fig.5. Bottom: 3D-shapes obtained integrating the normals[START_REF] Durou | Integration of a Normal Field in the Presence of Discontinuities[END_REF]. The normal fields n + and n -are only two solutions among 2 |Ω| , since any combination of both these normal fields is also plausible.

Figure 7 :

 7 Figure 7: The oversmoothing of the labeling obtained using high values of coefficients β c 2 , on the dataset from Fig. 6. From left to right: obtained labeling (black for n + , white for n -), XOR map between the estimated labels and the ground truth one (τ indicates the percentage of wrong labels), estimated normal field, and absolute angular error in degrees (MAE is the mean angular error). The choice β c 2 = β min indicates minimal coefficients (cf. (53)).

  83 • MAE = 13.65 • MAE = 11.58 • MAE = 0.43 • φ = 40 • MAE = 6.77 • MAE = 13.17 • MAE = 11.94 • MAE = 0.78 • φ = 60 • MAE = 12.64 • MAE = 15.87 • MAE = 15.76 • MAE = 2.06 •

Figure 8 :

 8 Figure8: Angular errors on the normals. The dataset from Fig.5was considered while increasing the zenithal angles φ for the lightings, in order to create stronger shadowing effects. Our method confines the errors in the shadow areas, while requiring neither a boundary condition nor parameter tuning (α 1 and α 2 are defined in Eq. (15)).

10 MAE = 1 .

 101 88 • MAE = 14.06 • MAE = 11.82 • MAE = 2.71 • d ρ = 50 MAE = 2.00 • MAE = 14.16 • MAE = 11.87 • MAE = 9.22 • d ρ = 100 MAE = 2.05 • MAE = 14.24 • MAE = 11.89 • MAE = 14.50 •

Figure 10 :

 10 Figure10: The 256 × 256 images were created with ρ everywhere equal to 1 except on a diagonal band with width d ρ pixels where it was set to 0.9. Assuming (wrongly) that ρ ≡ 1 induces a bias with our method, while existing methods are independent from the albedo.

Figure 12 :

 12 Figure12: Normals n + (x, y) and n -(x, y) move towards π as ρ(x, y) decreases from 1 towards a minimum value ρ inf (x, y) which corresponds to the limit case n + (x, y) = n -(x, y).

(a) I 1 (Figure 13 :Figure 14 :

 11314 Figure 13: (a-b-c) Three photographs of a plaster bust of Beethoven. (d) Albedo estimated using the PS3 technique, which is biased in the shadow areas. (e) Partition of Ω: Ω 3 (white), Ω 1,22 (red), and Ω 2,3 2 (blue). PS2 is applied in the Ω 2 sets, and PS3 elsewhere.

Figure 15 :

 15 Figure 15: Frontal (top) and side (bottom) views of the 3D-reconstructions obtained using PS3, and two versions of the combination PS3/PS2. Ours is able to recover fine-scale details more accurately in the shadow areas.

Figure 16 :

 16 Figure 16: First row: Three images of a Lambertian object, and the Ω 3 /Ω 2,3 2 partition. The area indicated in red lies almost entirely in the shadow in the first image. Second row: angular error (in degrees) obtained using, from left to right, the PS3 method (MAE = 11.84 • ), the PS2/PS3 combination from [9] (MAE = 11.27 • ), and the proposed one (MAE = 11.57 • ). The differential approach, which ensures smoothness, is globally more satisfactory than the non-differential one. Third row: 3D-reconstruction results. Fourth row: close-up on the estimated normals over the area indicated in red. The MAE are, respectively, 20.79 • , 25.49 • and 16.73 • . Our method achieves the best 3D-reconstruction in the shadow area.

Figure 17 : 2 ( 2 (

 1722 Figure 17: (a) An RGB image of a face illuminated by three directional, non-coplanar, color light sources. (b-c-d) Decomposition of the image (a) in three channels (red, green, blue). (e) Partition of Ω into sets Ω 3 (white), Ω R,B 2 (green), and Ω G,B 2 (blue).

5

  http://mi.eng.cam.ac.uk/research/projects/VideoNormals/

Figure 18 :

 18 Figure 18: 3D-reconstructions obtained from the single RGB image of Fig. 17-a. Only our method is able to restore the nostril area. The artifacts around the shadow area boundary are probably due to the overlapping spectra of the light sources: these artifacts are also visible in the images of Fig. 17-b and Fig. 17-c, although they should not appear.

  Figure 18: 3D-reconstructions obtained from the single RGB image of Fig. 17-a. Only our method is able to restore the nostril area. The artifacts around the shadow area boundary are probably due to the overlapping spectra of the light sources: these artifacts are also visible in the images of Fig. 17-b and Fig. 17-c, although they should not appear.

  Figure 18: 3D-reconstructions obtained from the single RGB image of Fig. 17-a. Only our method is able to restore the nostril area. The artifacts around the shadow area boundary are probably due to the overlapping spectra of the light sources: these artifacts are also visible in the images of Fig. 17-b and Fig. 17-c, although they should not appear.

The expressions[START_REF] Shi | A Benchmark Dataset and Evaluation for Non-Lambertian and Uncalibrated Photometric Stereo[END_REF] a posteriori explain the signification of the superscripts + and -.

It is stated in[START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF] that "the running time is nearly linear in practice".

http://www.ece.ncsu.edu/imaging/Archives/ImageDataBase/Industrial/


P Addesso

M Dalla Mura

L Condat

R Restaino

G Vivone

D Picone

J Chanussot

HYPERSPECTRAL PANSHARPENING USING CONVEX OPTIMIZATION AND COLLABORATIVE TOTAL VARIATION REGULARIZATION

Keywords: Data fusion, Hyperspectral pansharpening, Convex optimization, Total variation, Deconvolution

Hyperspectral pansharpening is a challenging research area and several methods have been recently developed to fuse low resolution hyperspectral and high resolution panchromatic images. In this paper we focus on a recent regularization method, called Collaborative Total Variation, exploiting a convex optimization algorithm. We evaluate the effectiveness of this novel approach in comparison to existing methods, and assess the performances on two datasets: a synthetic scene mimicking the characteristics of the Hyperion and ALI sensors and the Pavia University dataset.

INTRODUCTION

Images composed of several spectral channels allow to infer, with medium to low spatial resolution, the spectral signature of the materials present on a scene, as required by several applications, ranging from precision agriculture to mineralogy mapping [START_REF] Bioucas-Dias | Hyperspectral remote sensing data analysis and future challenges[END_REF]. On the other hand, high spatial resolution images can give a detailed representation of the surveyed area, but are in general less spectrally accurate (i.e., they feature broadband spectral responses). However, the availability of data characterized by both high spatial resolution and high spectral diversity is precluded in the real practice, due to physical constraints of acquisition devices. Hyperspectral Pansharpening [START_REF] Vivone | A critical comparison among pansharpening algorithms[END_REF][START_REF] Loncan | Hyperspectral pansharpening: A review[END_REF] addresses this limitation by generating a synthetic image through the fusion of a Hyper-Spectral (HS) image (with hundreds of bands in the Visible Near-Infrared and Shortwave Infrared) with a PANchromatic (PAN) image, characterized by a higher spatial resolution. Several algorithms have been proposed in the literature for performing hyperspectral pansharpening [START_REF] Loncan | Hyperspectral pansharpening: A review[END_REF]. They can be divided in three main families: Component Substitution (CS), MultiResolution Analysis (MRA) and Bayesian/Variational techniques. The two former groups are composed by classical pansharpening methods, extended to hyperspectral im-ages [START_REF] Vivone | A critical comparison among pansharpening algorithms[END_REF], whereas the latter can be seen as a particular instance of a deconvolution or deblurring problem.

Focusing on the latter family, variational techniques (e.g., based on convex optimization), have been profitably used in recent years [START_REF] He | A new pansharpening method based on spatial and spectral sparsity priors[END_REF][START_REF] Simões | A convex formulation for hyperspectral image superresolution via subspace-based regularization[END_REF]. In particular in [START_REF] Simões | A convex formulation for hyperspectral image superresolution via subspace-based regularization[END_REF] the authors have introduced the "HySure" method, in which an efficient implementation is attained by reducing the dimensionality of HS data and a regularization approach, based on a particular form of Vectorial Total Variation (VTV) [START_REF] Bresson | Fast dual minimization of the vectorial total variation norm and applications to color image processing[END_REF], is employed for dealing with the ill-posed nature of the problem. In this work we extend the approach proposed in [START_REF] Simões | A convex formulation for hyperspectral image superresolution via subspace-based regularization[END_REF] for the fusion of a HS and a PAN image by considering alternative regularization strategies based on the Collaborative Total Variation (CTV) paradigm. Different instances of CTV were defined and applied to denoising and deblurring problems in [START_REF] Duran | Collaborative total variation: A general framework for vectorial tv models[END_REF]. Therefore, in this paper we show that, according to the particular choice of the regularization term, the fused images exhibit different effects that directly impact the Hyperspectral Pansharpening results.

PANSHARPENING VIA TOTAL VARIATION

In this work we test several spatial regularization terms based on total variation, by relying on the framework proposed in [START_REF] Simões | A convex formulation for hyperspectral image superresolution via subspace-based regularization[END_REF]. To make the paper self-contained, we briefly describe it in the following, resorting to a vector notation. Namely, a multichannel image is represented by a matrix in which each row contains the lexicographically ordered pixels of a given band. More in detail, the HS image H ∈ R C×N HS , with C bands and N HS pixels, is assumed to be a spatially degraded version of the unknown high spatial and spectral resolution image Z ∈ R C×N P AN (with N P AN = ρ 2 N HS , where ρ > 1 is the resolution ratio), obtained according to the following model

H = ZBM + N h , (1)
where: i) B ∈ R N P AN ×N P AN is the spatial blurring matrix, representing the hyperspectral sensor's Point Spread Function (PSF), that is assumed to be band-independent; ii) M ∈ R N P AN ×N HS allows to perform a uniform subsampling of the Table 1: Synthetic dataset: performance obtained by the tested CTV norm averaged on 100 Monte Carlo trials. λ ϕ is the optimal value of the TV term weight, as defined in eq. ( 3) [START_REF] Simões | A convex formulation for hyperspectral image superresolution via subspace-based regularization[END_REF] 0.05 0.9919 1.8793 0.9334 0.9558 0.67 6.2949 11.438 0.7648 0.8799 39.9 image and gives rise to the lower spatial resolution of the hyperspectral image; iii) Z ∈ R C×N P AN is the high spatial and spectral resolution image (i.e., the objective of pansharpening); iv) N h ∈ R C×N HS is a matrix whose elements are samples of independent and identically distributed (i.i.d.) zeromean Gaussian noise with variance σ 2 HS . On the other hand, the PAN image P ∈ R 1×N P AN is supposed to be obtainable as:

SN R P AN = 40 dB, SN R HS = 30 dB SN R P AN = 20 dB, SN R HS = 20 dB Norm λ ϕ ERGAS SAM UIQI SCC λ ϕ ERGAS SAM UIQI SCC T [s] 2,2,1 (dbx)
P = RZ + N p (2)
where: i) R ∈ R 1×C is related to the Relative Spectral Response (RSR) of the Panchromatic band; ii) N p ∈ R 1×N P AN is a matrix whose elements are samples of independent and identically distributed (i.i.d.) zero-mean Gaussian noise with variance σ 2 P AN . Following [START_REF] Simões | A convex formulation for hyperspectral image superresolution via subspace-based regularization[END_REF], since Z lives in a subspace of dimensionality (significantly) lower than C, the models in (1) and ( 2) admit an alternative formalization. Specifically, Z is factorized as Z = EX, in which E is the set of basis (with cardinality lower than C) that spans the subspace of Z and X are the representation coefficients. This dimensionality reduction step, which can be obtained by using algorithms from the spectral unmixing literature [START_REF] Simões | A convex formulation for hyperspectral image superresolution via subspace-based regularization[END_REF], leads to a simplification of the problem that reflects in an increased robustness and eases the computation. Therefore, the aim is to estimate a pansharpened image Z (or equivalently its factorization EX) that gives the closest approximation of the HS and PAN im-age according to the models (1) and ( 2). This represents an ill-posed inverse problem that can be addressed by a convex optimization formalization, described by: minimize

X 1 2 H -EXBM 2 F + λm 2 P -REX 2 F + λ ϕ ϕ(X), (3)
where • F is the Frobenius norm and λ m = 1 [START_REF] Simões | A convex formulation for hyperspectral image superresolution via subspace-based regularization[END_REF]. Due to the complexity of the ill-posed problem, a regularization term ϕ(X) is required and its weight λ ϕ needs to be tuned. In the following we focus on this term, in order to identify the best choice according to the CTV theory.

COLLABORATIVE TOTAL VARIATION

The VTV regularization used in [START_REF] Simões | A convex formulation for hyperspectral image superresolution via subspace-based regularization[END_REF] is a particular instance of Total Variation on multivariate images (e.g., multi-and hyperspectral images). Indeed the authors of [START_REF] Duran | Collaborative total variation: A general framework for vectorial tv models[END_REF] showed that it is possible to build a framework, named CTV, in which the main idea is to choose the most appropriate regularization term ϕ(•) among the so-called Collaborative Norms (CNs), that operate on a multivariate image A = [A i,j,k ] ∈ R N ×C×M , where N is the number of the pixels, C is the number of the bands and M is the number of directional derivatives computed on each pixel.

More in detail, there are two general formulations for the CN's. The first one is based on the p,q,r norm. For example, if we associate the p norm to the derivative dimension (d), the q norm to the bands dimension (b) and the r norm to the pixel dimension (x), i.e. the application order is (dbx), we can obtain the p,q,r (dbx) norm :

||A|| p,q,r =    N i=1   C j=1 M k=1 |A i,j,k | p q/p   r/q    1/r . (4)
The other formulation is based on the Schatten p-norm (referred to as S p ) applied to two dimensions (usually spectra and derivatives) and then on the q norm for the third dimension (i.e. the pixel dimension), according to:

(S p , q ) (A) =    N i=1 A i,1,1 • • • A i,1,M . . . . . . . . . A i,C,1 • • • A i,C,M q S p    1/q , (5)
that defines the (S p (bd), q (x)) norm.

In our framework M = 2 and therefore we can put A in the form A = [XD h , XD v ], where XD h and XD v stand for the spatial derivatives of X in the horizontal and vertical directions, respectively. The range of possible choices is wide [START_REF] Duran | Collaborative total variation: A general framework for vectorial tv models[END_REF], due to the possibility of freely selecting the norm kind and the application order. In the following we will check the effectiveness of the most popular choices, such as the 2,2,1 (dbx) norm (used in "HySure" [START_REF] Simões | A convex formulation for hyperspectral image superresolution via subspace-based regularization[END_REF]) and the (S 1 (bd), 1 (x)) norm (the well-known nuclear norm).

EXPERIMENTAL RESULTS

The effect of selecting different CNs is evaluated via a reduced resolution assessment [START_REF] Wald | Fusion of satellite images of different spatial resolutions: Assessing the quality of resulting images[END_REF] on two different datasets: a synthetic one, designed to mimic the acquisitions of the Hyperion and ALI sensors, and the well-known Pavia University dataset. The fusion results are evaluated via several in- dices [START_REF] Vivone | A critical comparison among pansharpening algorithms[END_REF]9,[START_REF] Zhou | A wavelet transform method to merge landsat TM and SPOT panchromatic data[END_REF], i.e.: i) The Erreur Relative Globale Adimensionelle de Synthèse (ERGAS) that is a generalization of the Root Mean Square Error (RMSE); ii) The Spectral Angle Mapper (SAM), useful to evaluate the spectral distortion; iii) The Universal Image Quality Index (UIQI), that takes into account the correlation, the error and the contrast between the fused image and the reference image; iv) The Spatial Correlation Coefficient (SCC) computed between the details of the reference image and the fused one, aimed to assess the correct rendering of the details in the fused image.

Synthetic dataset

This dataset has been created by using a collection of ten geometric shapes: each of them is associated to a spectral signature of a different material. The background, divided in four quarters, is composed by a mixing (different for each quarter) of the previous materials. Three images have been generated (see Fig. 1 (a)-(c)): i) A HS image, playing the role of the Ground Truth (GT), obtained by using the Relative Spectral Response (RSR) of the Hyperion sensor; ii) A PAN image, obtained using the RSR of the ALI sensor. Moreover, white Gaussian noise is added in order to have a Signal to Noise Ratio (SNR) equal to a fixed value SN R P AN ; iii) A HS image, Table 2: Pavia University dataset: performance obtained by some CTV, CS and MRA algorithms, averaged on 100 Monte Carlo trials. λ ϕ is the optimal value of the TV term weight in eq. ( 3). 1, two different norms achieve the best results in the high SNR and low SNR scenarios. The 2,1,1 (bdx) norm represents the best option in the first case (SN R P AN = 40 dB, SN R HS = 30 dB) and the nuclear norm is the optimal choice (S 1 (bd), 1 (x)) in the second setting (SN R P AN = 20 dB, SN R HS = 20 dB). A significant improvement is achieved w.r.t. the 2,2,1 (dbx) norm used in "HySure" [START_REF] Simões | A convex formulation for hyperspectral image superresolution via subspace-based regularization[END_REF] in terms of SAM and ERGAS. Focusing on the spatial distribution of the SAM, reported in Fig. 2, we can see that the improvements are mainly related to the borders of the geometrical shapes, that are sharper for 2,1,1 (bdx) and (S 1 (bd), 1 (x)) (see, for instance, the little square in the bottom-left corner). Noticeably, CTV algorithms based on these three norms and on 1,1,1 (bdx) (that obtains appreciable SAM values in the high noise scenario) have comparable (and relatively low) computation times T. The use of other norms results in higher computation times and lower performance, thus constituting worse choices.

SN R

Pavia University dataset

This dataset is based on a widely used HS image acquired by the ROSIS sensor. Also in this case, three images were created (see Fig. 3 (a)-(c)): i) The original ROSIS image plays the role of GT; ii) The PAN image is simulated by using the RSR of IKONOS sensor; iii) The HS image is obtained by blurring the GT via a Starck-Murtagh filter [START_REF] Strang | Wavelets and Filter Banks[END_REF] and downsampling by a factor ρ = 4. Similarly to the other dataset, both PAN and HS images are corrupted by adding white Gaussian noise in order to obtain the desired values of the SN R P AN and SN R HS , respectively. In this case, we considered for comparison only the fastest CTV-based algorithms, i.e., the ones based on 2,2,1 (dbx), 1,1,1 (bdx),

2,1,1 (bdx) and (S 1 (bd), 1 (x)) norms. Moreover, for sake of comparison, we assess the performance of some conventional pansharpening techniques belonging to both the CS class (i.e., GS and GSA, based on the Gram-Schmidt procedure, and the Principal Component Analysis (PCA) method) and to the MRA class (i.e., the High-Pass Filtering (HPF) and the A-Trous Wavelet Transform (ATWT) method) [START_REF] Vivone | A critical comparison among pansharpening algorithms[END_REF]. Looking at the results shown in Table 2, we see that CTV-based algorithms outperform the other techniques. More specifically, we can see that 2,1,1 (bdx) yields a poor performance: indeed, abrupt variations (caused for example by the shadows) are present in all the bands of this image. The 2,1,1 (bdx) norm is not able to reduce spectral artifacts on the corresponding edges [START_REF] Duran | Collaborative total variation: A general framework for vectorial tv models[END_REF][START_REF] Duran | On the Implementation of Collaborative TV Regularization: Application to Cartoon+Texture Decomposition[END_REF]. On the other hand, the nuclear norm (S 1 (bd), 1 (x)) still outperforms the "HySure" 2,2,1 (dbx) one, in particular in terms of SAM, both for low and high noise. Finally, as in the synthetic dataset, 1,1,1 (bdx) norm shows good performance in the high noise scenario.

CONCLUSIONS

In this work we compared the effectiveness of several norms to implement the Collaborative Total Variation framework for the fusion of HS and PAN images. The numerical results obtained with two different datasets indicate that the nuclear norm S 1 (bd), 1 (x)) allows to obtain superior performance with respect to the 2,2,1 (dbx) norm used in the original "HySure" algorithm. Future work will be devoted to the application of these algorithms to other pansharpening scenarios, e.g., the fusion of HS and MS images and of MS and PAN ones.

Fig. 1 :

 1 Fig. 1: Synthetic dataset: GT image (a), PAN image (b), HS image(c), example of fused image using 2,2,1 (dbx) norm (d).

Fig. 2 :

 2 Fig. 2: Synthetic dataset: SAM map for 2,2,1 (dbx) norm (a), 1,1,1 (bdx) norm (b), 2,1,1 (bdx) norm (c) and (S 1 (bd), 1 (x)) norm (d). Relevant parameters are SN R P AN = 40 dB and SN R HS = 30 dB.

Fig. 3 :

 3 Pavia University dataset: GT image (a), PAN image (b), HS image(c), Fused image (using 2,2,1 (dbx) norm) (d).

  P AN = 40 dB, SN R HS = 30 dB SN R P AN = 20 dB, SN R HS = 20 dB degrading the GT according to the PSF of the Hyperion sensor: the resolution ratio is ρ = 3. Also in this case, white Gaussian noise is added in order to have an SNR equal to a fixed value SN R HS .As shown by the numerical results reported in Table

	Algorithm	λ ϕ	ERGAS SAM	UIQI	SCC	λ ϕ	ERGAS SAM	UIQI	SCC
	EXP	-	7.3828	5.2903 0.7686 0.5517	-	7.4267	5.7044 0.7623 0.5533
	HPF	-	5.8474	7.0084 0.8792 0.7261	-	6.8810	8.0241 0.8105 0.7345
	ATWT	-	5.9789	8.0134 0.8789 0.7404	-	6.9640	8.8917 0.8154 0.7469
	GS	-	5.3330	6.3522 0.8839 0.7411	-	6.1975	7.0914 0.8276 0.7433
	GSA	-	6.0100	9.2121 0.8798 0.7409	-	7.4846	10.501 0.8056 0.7460
	PCA	-	7.4006	9.3170 0.7843 0.6854	-	8.0904	9.9096 0.7387 0.6899
	CTV: 2,2,1 (dbx) [5]	0.002	3.8160	4.8204 0.9411 0.7804	0.02	4.2458	5.0985 0.9168 0.7494
	CTV: 1,1,1 (bdx)	0.002	3.8431	4.8398 0.9387 0.7777 0.005	4.2891	5.0138 0.9147 0.7679
	CTV: 2,1,1 (bdx)	0.002	3.9325	4.9662 0.9370 0.7763	0.01	4.3937	5.2126 0.9116 0.7600
	CTV: (S 1 (bd), 1 (x)) 0.002	3.7809	4.7396 0.9421 0.7801	0.02	4.2097	5.0285 0.9179 0.7436
	obtained by spatially


iMacs deficiency in nod mice : implication in cyclophosphamide_induced diabetes ? Stéphanie Ferret-Bernard, Pierre Saï, Jean-Marie Bach To cite this version: Stéphanie Ferret-Bernard, Pierre Saï, Jean-Marie Bach. iMacs deficiency in nod mice : implication in cyclophosphamide_induced diabetes ?. Congrès Annuel du Club Francophone des Cellules Dendritiques (CFCD), Dec 2002, Paris, France. 2002, Congrès 2002. hal-01399474