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

François Lauze
email: francois@diku.dk

Jean-Denis Durou

A L 1 -TV Algorithm for Robust Perspective Photometric Stereo with Spatially-Varying Lightings

Keywords: Uncalibrated photometric stereo, Spatially-varying lightings, Perspective projection, Total variation, Proximal methods

We tackle the problem of perspective 3D-reconstruction of Lambertian surfaces through photometric stereo, in the presence of outliers to Lambert's law, depth discontinuities, and unknown spatiallyvarying lightings. To this purpose, we introduce a robust L 1 -TV variational formulation of the recovery problem where the shape itself is the main unknown, which naturally enforces integrability and permits to avoid integrating the normal field.

Introduction

Photometric stereo refers to the problem of inferring the shape of an object, given a set of m images of this object acquired under different known (calibrated photometric stereo) or unknown (uncalibrated photometric stereo) illuminations but from the same point of view, by inverting a generative image model accounting for the projection model (orthographic or perspective), the object reflectance (Lambert, Phong, etc.) and the lighting model (directional, punctual, etc.). Early photometric stereo considered orthographic projection, Lambertian reflectance and directional lightings [START_REF] Woodham | Photometric method for determining surface orientation from multiple images[END_REF], which is the simplest case. Extending photometric stereo to more general models has become an important research direction, along with the natural need for fast and robust solutions. For example, perspective projection can replace the orthographic projection in both calibrated [START_REF] Mecca | A direct differential approach to photometric stereo with perspective wiewing[END_REF] and uncalibrated [START_REF] Papadhimitri | A new perspective on uncalibrated photometric stereo[END_REF] cases, the Lambertian assumption can be relaxed so as to take specular highlights into account [START_REF] Georghiades | Incorporating the Torrance and Sparrow model of reflectance in uncalibrated photometric stereo[END_REF] and robustness to outliers like shadows can be obtained through numerous techniques such as sparse regression [START_REF] Ikehata | Robust photometric stereo using sparse regression[END_REF]. However, extending photometric stereo to non-directional lightings remains an open and challenging problem. We focus on this aspect of the problem and propose a robust variational formulation of perspective photometric stereo with spatially-varying lightings, considered as additional unknowns so as not to impose a parametric lighting model. The different hypotheses our approach relies on are described hereafter.

Image Formation Model

Radiance Model. Lambert's law assumes that the surface S reflects light diffusively, the reflectance at x ∈ S being characterised by the albedo ρ(x) ∈ [0, 1]. Given a parallel uniform light beam s 0 ∈ R 3 , oriented towards the source, and the unit outward normal n(x) ∈ R 3 to the surface at x, the emitted Lambertian radiance r e (x) at x is proportional to the irradiance max {0, s 0 • n(x)}, where the max operator aims at modeling the self-shadows. Forgetting this operator (shadows will be dealt with by robust estimators), the Lambertian emitted radiance reads:

r e (x) = ρ(x) π s 0 • n(x) (1)
In this work, we do not suppose that the lightings consist in parallel uniform beams. Even more, we do not use any explicit lighting model, allowing us to cope with, for instance, near point light sources, large sources, self-and cast-shadows and even secondary reflections. In the most general case, the irradiance at x is equal to the sum of elementary contributions coming from all the unit directions v(θ, φ) characterised by the spherical angles (θ, φ) defined with respect to n(x):

r e (x) = ρ(x) π θ∈[0,π/2], φ∈[0,2π] r i (x, θ, φ)v(θ, φ) • n(x) sin θ dθ dφ (2)
where r i (x, θ, φ) is the incident radiance at x in the direction v(θ, φ). Let s(x) = θ,φ r i (x, θ, φ)v(θ, φ) sin θ dθ dφ. We call the vector field s the light field. Eq. ( 2) then provides a more general Lambertian emitted radiance model than [START_REF] Alldrin | Resolving the generalized bas-relief ambiguity by entropy minimization[END_REF]:

r e (x) = ρ(x) π s(x) • n(x) (3)
In the photometric stereo framework, m different light fields s i , i ∈ [1, m], successively illuminate the scene, providing m images I i , i ∈ [1, m], assumed to be graylevel images. The graylevel is assumed to be proportional to the emitted radiance r e (x) as expressed in (3). To our knowledge, this very general radiance model has not been used in the context of photometric stereo so far. In most applications, a single point light source at infinity is assumed, so that each lighting is represented by a constant vector in R 3 [START_REF] Woodham | Photometric method for determining surface orientation from multiple images[END_REF][START_REF] Georghiades | Incorporating the Torrance and Sparrow model of reflectance in uncalibrated photometric stereo[END_REF][START_REF] Alldrin | Resolving the generalized bas-relief ambiguity by entropy minimization[END_REF][START_REF] Papadhimitri | A new perspective on uncalibrated photometric stereo[END_REF][START_REF] Quéau | Solving uncalibrated photometric stereo using total variation[END_REF]. Closer to this general expression is the important work of Basri et al. in [START_REF] Basri | Photometric stereo with general, unknown lighting[END_REF], where a distribution of point light sources at infinity is assumed, allowing the authors to develop shape recovery algorithms based on low order spherical harmonic decompositions.

Camera Model. We assume the camera is a calibrated pinhole camera, represented by its focal length f and its principal point, considered as the origin of the image domain Ω ⊂ R 2 . Each pixel x p ∈ Ω is then uniquely associated to a 3D point x(x p ) = [x(x p ), y(x p ), z(x p )] ⊤ ∈ R 3 , where z > 0 is called the depth, so that the surface of the object to be reconstructed can be represented as the set of 3D points S = {x(x p ), x p ∈ Ω}, with:

x(x p ) = exp(u(x p )/f ) 1 f x p 1 , u(x p ) = f log z(x p ) (4)
It follows from (4) that the unit outward normal n(x) to the surface at x can be equivalently defined on the image domain Ω, through (see e.g. [START_REF] Mecca | A direct differential approach to photometric stereo with perspective wiewing[END_REF]):

n u (x p ) = 1 ∇u (x p ) 2 + 1 + 1 f x p • ∇u(x p ) 2 -∇u (x p ) 1 + 1 f x p • ∇u (x p ) (5)
In the radiance model (3), the albedo ρ and the light field s are defined on the surface S, which is a subset of R 3 . Yet, since they are independent from the projection model, they can be equivalently defined onto the image domain Ω through the mapping (4). We will thus refer to their values as ρ(x p ) and s(x p ). Forgetting the normalisation factor 1 π , the image formation model finally reads:

I i (x p ) = ρ(x p ) s i (x p ) • n u (x p ), x p ∈ Ω, i ∈ [1, m] (6)
knowing that this equality is in fact a relation of proportionality.

Overview of our Contributions

In most of literature on photometric stereo, the main unknown is the couple (albedo, normal field) which is estimated from the set of linear equations ( 6) without any prior on the camera model: knowledge of the camera parameters only has importance when integrating this estimated normal field into a depth map [START_REF] Durou | Integrating the normal field of a surface in the presence of discontinuities[END_REF]. Yet, recent photometric approaches explicitely benefit from the perspective camera model: it is shown in [START_REF] Papadhimitri | A new perspective on uncalibrated photometric stereo[END_REF] that the directional uncalibrated photometric stereo problem is better constrained under perspective projection, provided the normal field is considered integrable. Also, Mecca et al. show in [START_REF] Mecca | A direct differential approach to photometric stereo with perspective wiewing[END_REF] that the system of nonlinear PDEs in z = exp(u/f ) resulting from ( 5) and ( 6) can be solved through a semi-Lagrangian scheme, avoiding the classical two-steps approach. The path we follow here combines the benefits from both these recent works, and improve their robustness through the introduction of a variational formalism.

To emphasise the advantages of our approach over prior work, let us state that: 1) it is the first variational approach to uncalibrated photometric stereo; 2) it enhances robustness through L 1 -TV optimisation (only least-squares approaches were considered in the uncalibrated case so far [START_REF] Georghiades | Incorporating the Torrance and Sparrow model of reflectance in uncalibrated photometric stereo[END_REF][START_REF] Papadhimitri | A new perspective on uncalibrated photometric stereo[END_REF][START_REF] Quéau | Solving uncalibrated photometric stereo using total variation[END_REF]); 3) depth discontinuities can be recovered without the need for a posteriori robust normal field integration through e.g., techniques like those described in [START_REF] Durou | Integrating the normal field of a surface in the presence of discontinuities[END_REF]; and 4) unknown, spatially-varying, lightings and albedo are considered (the one-step approach from [START_REF] Mecca | A direct differential approach to photometric stereo with perspective wiewing[END_REF] assumes calibrated directional lightings and known albedo).

The rest of this paper is organised as follows. Using a Bayesian rationale, we use the hypotheses above to derive in Section 2 a generic maximum a posteriori formulation of the recovery problem, and a variational formulation where regularisations use total variation semi-norms. Then, we propose in Section 3 a proximal algorithm for minimising the associated energy, before demonstrating in Section 4 the benefit of using spatially-varying lightings and the proposed robust recovery.

We derive in this section a variational approach for inverting the image formation model ( 6), via ideas from Bayesian inference for factoring the posterior probability, in the spirit of Mumford's work [START_REF] Mumford | Bayesian rationale for the variational formulation[END_REF]. For the sake of compactness, we will denote the set of images i = [I 1 , . . . , I m ] ⊤ (vector field Ω → R m ) and that of lightings S = s 1 , . . . , s m (matrix field Ω → R 3×m ).

MAP Estimation

We first propose to recover (u, ρ, S) as the maximum a posteriori (MAP) of the distribution

P(u, ρ, S|i) = P(i|u, ρ, S) P(u, ρ, S) P(i) (7)
where P(i) is the evidence, which can be discarded since it is constant and plays no role in MAP, P(i|u, ρ, S) is the likelihood, and P(u, ρ, S) is the prior, which is factored as P(u, ρ, S) = P(S|u, ρ)P(u, ρ). In all photometric stereo works, u and ρ are expressed as independent functions on Ω, which is also the point of view used here, thus we factor P(u, ρ) = P(u)P(ρ).

For sake of simplicity, we assume that the lightings are independent from the surface and the albedo as well, so that the conditional prior simplifies: P(S|u, ρ) = P(S). Of course, this simplification is abusive in the presence of secundary reflections. Yet, such reflections are in general sufficiently sparse to be neglected almost everywhere, and their effect on depth recovery can be limited by an appropriate choice of regularisation, as the one proposed hereafter. Studying the benefit of using a joint conditional prior is left as a future research direction.

With these assumptions, we eventually obtain P(u, ρ, S) = P(u)P(ρ)P(S). Maximisation of the probability ( 7) is then equivalent to minimisation of the neg-log-posterior E(u, ρ, S) = -log P(u, ρ, S|i) -log P(i) which is given by: E(u, ρ, S) = -log P(i|u, ρ, S) -log P(u) -log P(ρ) -log P(S) [START_REF] Goldluecke | The natural vectorial total variation which arises from geometric measure theory[END_REF] where the first term is a data term, and the others are regularisation measures.

Continuous Variational Energy

Data Term. Lambert's law (6) does generally not hold perfectly: noise in the measurements will prevent the equality. We assume that Lambert's law residuals are iid Laplacian, with some outliers that will account for the shadows (unexplained low graylevels), the highlights (unexplained high graylevels) and the depth discontinuities (∇u not defined). We thus use a L 1 -norm based data term that accounts for the sparsity of such effects:

E D (u, ρ, S)= m i=1 Ω I i (x p )-ρ(x p )s i (x p ) •n u (x p ) dx p = m i=1 I i -ρs i • n u L 1 (Ω) (9)
Regularisations. In the 3D-reconstruction framework, surfaces are usually assumed to be differentiable, but this hypothesis is rarely realistic because of the presence of edges and depth discontinuities. At a macroscopic scale, surfaces should, more reasonably, be supposed to be smooth almost everywhere, with continuous, non-differentiable edges and some depth discontinuities. A regulariser that naturally allows for the corresponding class of depth functions is the total variation (TV) semi-norm:

J(u) = Ω ∇u(x p ) dx p = ∇u L 1 (Ω) (10)
As suggested in literature [START_REF] Alldrin | Resolving the generalized bas-relief ambiguity by entropy minimization[END_REF][START_REF] Quéau | Solving uncalibrated photometric stereo using total variation[END_REF], we restrict our study to the case of piecewisesmooth ("Arlequin-like") albedos, which can be enforced by a TV semi-norm as well:

J(ρ) = Ω ∇ρ(x p ) dx p = ∇ρ L 1 (Ω) (11)
In the directional orthographic case, the problem is inherently ill-posed, but it is shown in [START_REF] Quéau | Solving uncalibrated photometric stereo using total variation[END_REF] that TV-regularisation of both the depth and the albedo reduces the ambiguities to a simple translational ambiguity. The perspective case is better constrained (u can be recovered up to a scale ambiguity on S without regularisation [START_REF] Papadhimitri | A new perspective on uncalibrated photometric stereo[END_REF]), so these regularisations would essentially enforce smoothness. Well-posedness is less clear in perspective case with spatially-varying lightings. Yet, it seems rather intuitive that, by limiting the variations of the lightings, the problem gets closer to the directional one, at least piecewise. This general discussion is left as future work, but it invites us to introduce a regularisation measure for the lightings as well. Since the depth discontinuities of both u and s i will, in general, coincide, it is reasonable to assume that s i should have bounded variations too, and to introduce its TV-regularisation:

J(s i ) = Ω J(s i )(x p ) F dx p = J(s i ) L 1 (Ω) , i ∈ [1, m] (12)
where J(s i ) is the Jacobian matrix of s i and • F denotes the Frobenius norm.

For the sake of compactness, we will denote J(S) = m i=1 J(s i ). As stated earlier, the independence of S and u serves as a simplifying hypothesis, but it has no physical motivation: in future work, a joint regularisation of both these fields shall be introduced, using for instance a coupled L 1 semi-norm [START_REF] Goldluecke | The natural vectorial total variation which arises from geometric measure theory[END_REF].

Energy. A continuous analogue of the MAP problem of minimising the neg-log posterior (8) can thus be written as the minimisation of E(u, ρ, S) = E D (u, ρ, S)+ αJ(u) + βJ(ρ) + γJ(S), with α, β, γ some positive weights. This functional is however not coercive in u, because of the scale ambiguity on S inherent to monocular perspective 3D-reconstruction. Coercivity can be forced by setting arbitrarily the depth in one point, or by adding a quadratic prior u 0 :

E(u, ρ, S) = m i=1 I i -ρ s i • n u L 1 + ǫ 2 u -u 0 2 L 2 +αJ(u)+βJ(ρ)+γJ(S) (13)
with ǫ > 0 (u 0 can for instance be a prior on the mean camera-object distance).

Proximal Recovery Using Split-Bregman Iterations

We propose a proximal algorithm [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] for dealing with the TV semi-norms in [START_REF] Osher | An iterative regularization method for total variation-based image restoration[END_REF], and introduce Split-Bregman iterations for decoupling the L 1 data term.

Adopting the same notations as in [START_REF] Goldstein | The Split Bregman method for L1-regularized problems[END_REF] for the Bregman variables (d, b), whose components will be denoted (d i , b i ), the augmented energy is then defined as:

Ē(u, ū, ρ, ρ, S, S, d, b) = m i=1 d i L 1 + 1 2θ d m i=1 d i -I i -ρ si •n ū -b i 2 L 2 + ǫ 2 ū-u 0 2 L 2 + 1 2θ u ū-u 2 L 2 +αJ(u)+ 1 2θ ρ ρ-ρ 2 L 2 +βJ(ρ)+ 1 2θ S S-S 2 L 2 + γJ(S) (14)
for some small positive θ d , θ u , θ ρ , θ S . We minimise the augmented energy ( 14) iteratively w.r.t. each variable, according to Algorithm 1.

Algorithm 1 L 1 -TV Perspective Photometric Stereo

Input:

Images I i , i ∈ [1,
; ρ0 = ρ 0 , ū0 = u 0 ; d i,0 = b i,0 = 0, si,0 = s i,0 , i ∈ [1, m]. for k = 1 . . . N do ūk+1 = argmin ū 1 2θ d m i=1 d i,k -I i -ρ k si,k •nū -b i,k 2 L 2 + ǫ 2 ū-u0 2 L 2 + 1 2θu ū-u k 2 L 2 (15)
ρk+1 = argmin ρ 1 2θ d m i=1 d i,k -I i -ρ si,k •n ūk+1 -b i,k 2 L 2 + 1 2θρ ρ-ρ k 2 L 2 (16)
si,k+1 = argmin

si 1 2θ d d i,k -I i -ρ k+1 si •n ūk+1 -b i,k 2 L 2 + 1 2θ S si -s i,k 2 L 2 , i ∈ [1, m] (17)
d i,k+1 = argmin d i d i L 1 + 1 2θ d d i -I i -ρ k+1 si,k+1 •n ūk+1 -b i,k 2 L 2 , i ∈ [1, m] (18) b i,k+1 = b i,k + I i -ρ k+1 si,k+1 •n ūk+1 -d i,k+1 , i ∈ [1, m] (19)
u k+1 = argmin u 1 2θu ūk+1 -u 2 L 2 + αJ(u) (20)
ρ k+1 = argmin ρ 1 2θρ ρk+1 -ρ 2 L 2 + βJ(ρ) (21)
s i,k+1 = argmin s i 1 2θs si,k+1 -s i 2 L 2 + γJ(s i ), i ∈ [1, m] (22)
end for x = exp(u/f ) xp ⊤ /f, 1 ⊤

Solutions of Convex Subproblems in Algorithm 1

The update (18) in d is a basis pursuit problem, which is solvable explicitely by shrinkage:

d i,k+1 = shrink(I i -ρk+1 si,k+1 • n ūk+1 + b i,k , 1/θ d ), i ∈ [1, m] (23)
with shrink(x, γ) = sign(x) max{|x|-γ, 0}. The update (19) in b is the Bregman update. Justification for this update can be found in [START_REF] Goldstein | The Split Bregman method for L1-regularized problems[END_REF]. It is the equivalent of "adding back the noise" in iterative image denoising methods [START_REF] Osher | An iterative regularization method for total variation-based image restoration[END_REF].

Regarding the updates in u, ρ and S (Eqs (20), (21), and ( 22)), they are instances of the L 2 -TV problem which we solve using Chambolle's dual projection algorithm [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF], with the natural vectorial extension presented in [START_REF] Bresson | Fast dual minimization of the vectorial total variation norm and applications to color image processing[END_REF] in the case of S. We experimentally noticed that the convergence rate is better when only performing a few projection iterations (typically 5) at each global iteration k, rather than solving these problems to full convergence (the same behaviour is observed in [START_REF] Goldstein | The Split Bregman method for L1-regularized problems[END_REF]). As for the update ( 16) in ρ, writing the normal equations leads to the explicit update:

ρk+1 = θ d ρ k + θ ρ m i=1 (I i -d i,k + b i,k )(s i,k • n ūk+1 ) θ d + θ ρ m i=1 (s i,k • n ūk+1 ) 2 (24)
In a similar way, from the vectorial normal equations derived from [START_REF] Woodham | Photometric method for determining surface orientation from multiple images[END_REF], one obtains the explicit solution in si , i ∈ [1, m]:

si,k+1 = θ d I 3 +θ S ρk+1 2 n ūk+1 n ⊤ ūk+1 -1 θ d s i,k +θ S ρk+1 (I i -d i,k +b i,k )n ūk+1 (25)
which can be computed e.g., through Cholesky factorisation.

Solution in ū

We now describe the minimisation in ū with greater care. Let us denote:

E k ū (ū)= 1 2θ d m i=1 d i,k I i -ρ k si,k •n ū -b i,k 2 L 2 + ǫ 2 ū-u 0 2 L 2+ 1 2θ u ū-u k 2 L 2 (26)
Proposition 1. The necessary optimality condition for ū reads as the vanishing of the gradient of the energy (26) w.r.t. ū, which is given by:

∇ ūE k ū (ū) =div f k ū θ d (∇ū + (1 + ∇ū • u)u) + div 1 θ d π 1,2 (g k ū) -div π 3 g k ū θ d u + ǫ + 1 θ u ū -ǫu 0 - 1 θ u u k (27)
where π 1,2 and π 3 are the projections

π 1,2 (a, b, c) = (a, b), π 3 (a, b, c) = c, and
f k ū = m i=1 Ĩi,k Ĩi,k +d i,k -I i -b i,k ∇ū 2 + (1 + ∇ū • u) 2 , g k ū = m i=1 ρk Ĩi,k +d i,k -I i -b i,k ∇ū 2 + (1 + ∇ū • u) 2 si,k u = x p /f, Ĩi,k = ρk si,k • n ū (28)
Proof. The directional derivative of n ū at ū in the direction w is given by:

d w n ū(ū) = - ∇w• (∇ū+(1+∇ū• u)u) ∇ū 2 +(1+∇ū• u) 2 n ū- 1 ∇ū 2 +(1+∇ū• u) 2 ∇w -∇w• u (29) Let D i (ū) = 1 2θ d d i,k -I i -ρ k si,k •n ū -b i,k 2 L 2 .
Using the chain rule, its first variation L i (w) is given by:

L i (w) = - 1 θ d Ω f i,k ū ∇w • (∇ū + (1 + ∇ū • u)u) - 1 θ d Ω g i,k ū • ∇w -∇w • u ( 30
)
where f i,k ū and g i,k ū represent the terms inside the sums in (28). The first variation of ( 26) is thus given by:

L(w) = - 1 θ d Ω m i=1 f i,k ū (∇ū+(1+∇ū • u)u)•∇w+ Ω ǫ+ 1 θ u ū-ǫu 0 - 1 θ u u k w - 1 θ d Ω π 1,2 m i=1 g i,k ū -π 3 m i=1 g i,k ū u • ∇w (31)
Using Dirichlet boundary conditions on ū, application of the Green formula eventually provides the result announced in (27).

⊓ ⊔

This ressembles the steady state of a reaction-diffusion equation. The data term of (26) being nonlinear and non-convex in ū, there is no guarantee for the negativity of f k ū . For this reason, we choose a descent solver with a quasi-Newton step (BFGS) in the implementation. We now evaluate experimentally the benefit of using the proposed approach.

Experiments

Robustness to Noise and Outliers. We first evaluate the ability of the proposed proximal algorithm to enhance the results of classical photometric stereo. We create 20 images, of size 128 × 128, of a "Canadian tent" shape with a "pears" albedo, illuminated from m = 20 known lightings, and corrupt them simultaneously by an additive zero-mean Gaussian noise with standard deviation equal to σ% of the maximum graylevel, and a salt-and-pepper noise affecting p% of the pixels. We use the method from [START_REF] Ikehata | Robust photometric stereo using sparse regression[END_REF] to estimate the initial albedo ρ 0 and the normals which minimise the L 1 data term [START_REF] Goldstein | The Split Bregman method for L1-regularized problems[END_REF]. Up to this point, in view of the piecewise smooth nature of the shape, robust integration should be considered to obtain a depth map [START_REF] Durou | Integrating the normal field of a surface in the presence of discontinuities[END_REF]. To emphasise the ability of our scheme to recover such shapes as well, we integrate the normals using the DCT solver from [START_REF] Simchony | Direct analytical methods for solving Poisson equations in computer vision problems[END_REF] instead, with homogeneous Dirichlet BC, before iteratively refining the depth and the albedo through the proposed L 1 -TV scheme (Figure 1): this allows recovering the depth discontinuities, while simultaneously denoising the albedo. We used θ d = 0.1, θ u = 0.1θ d , θ ρ = 0.001θ u , θ S = θ ρ , ǫ = 10 -8 , u 0 = 0, α = 2/θ u , β = 5α, and γ = α. [START_REF] Ikehata | Robust photometric stereo using sparse regression[END_REF], using a least squares integrator with homogeneous Dirichlet BC; Fifth column: results after L 1 -TV refinement. State-ofthe-art methods provide a rather noisy albedo, and the quality of the final depth map is widely dependent from the choice of a specific integrator. Rather than choosing a robust integrator, as advised e.g. in [START_REF] Durou | Integrating the normal field of a surface in the presence of discontinuities[END_REF], we believe that it makes more sense to use the proposed scheme to directly refine the depth map and the albedo from the images, so as to avoid the propagation of any bias due to normal estimation. Indeed, rather than seeking a piecewise-smooth shape explaining a possibly biased normal field, it makes more sense to look for a piecewise-smooth shape explaining the images themselves.

Influence of the parameters. Since the proposed proximal algorithm involves numerous parameters, it is natural to question their influence. Let us first make a distinction between the model parameters α, β, γ and ǫ, and the optimisation parameters θ ⋆ . Those latter will mainly affect the convergence rate (though, since we consider nonconvex optimisation, they may also affect the results), while the choice of the former corresponds to some a priori knowledge about the regularity of the unknowns. If we choose high values for α, β, γ, we will get oversmoothed results, while the 3D-reconstruction may not be robust to noise, or not well-conditioned, if we choose low values. The parameter ǫ only enforces coercivity, thus any low value should be satisfactory (we set ǫ = 10 -8 in all the experiments). As for the optimisation parameters, setting high θ ⋆ values will result in slow convergence, while the optimisation may become unstable with low values. In our tests, we experimented no major difficulty in finding a set of these numerous parameters offering "reasonable" results. Tuning more precisely one or the other parameter can be a little more tricky, as illustrated in Figure 2. We used the same data as in the previous test, with σ = 0.1 and p = 0.05, for studying the evolution of the RMSE in u (the scale ambiguity on S is a posteriori solved to minimise this RMSE) with respect to a specific choice of α, θ d or θ u (the other parameters behave similarly). Finding heuristics for automatic selection of the parameters, based on these observations, is left as future work. Experimental Setup and Lighting Initialisation. For tests on real-world data, we used a calibrated camera surrounded by twenty LEDs oriented towards a white-painted two cents coin. The LEDs are only controllable by groups of 5, providing 15 different lighting configurations with either 5, 10, 15 or 20 LEDs turned on together. Thanks to the proposed spatially-varying lightings, we do not need model-driven calibration of the lightings: to initialise the light fields, we used a planar calibration grid composed of regularly located white truncated pyramids, before interpolating them so as to obtain dense light fields. We first evaluate the benefit of using spatially-varying lightings, by comparing the results of the classical two-steps approach, using either directional lightings or the calibrated spatially-varying ones (least-squares estimations are used). Results shown in Figure 3 prove that the widely-spread idea according to which photometric stereo can recover high-frequency details, but not lowfrequency ones is wrong: it is only a matter of correctly estimating the lightings.

L 1 -TV Refinement from Rough Initialisation. We now question the importance of the initialisation for the full (shape, reflectance and lightings) proximal recovery, regarding the non-convexity of the energy. Due to the computational cost required to perform one Cholesky factorisation per pixel at each iteration in the lightings estimation, we consider a small close-up on the "E" letter of the coin, with size 201 × 146 pixels. The evolution of the anisotropic diffusion process, starting from a flat shape with uniform albedo and directional lightings, is illustrated in Figure 4: it proves that the proposed proximal algorithm converges towards a realistic solution even if a rough initialisation is considered. We used θ d = 0.05, θ u = 100θ d , θ ρ = θ d, θ S = 0.1θ d , ǫ = 10 -8 , α = 0.05/θ u , β = 0.05/θ ρ , u 0 = 0 and γ = 0.0001/θ S . We proposed a variational formulation for the joint recovery of shape, reflectance and spatially-varying lightings from observed images, derived from a formal MAP approach. Despite the apparent difficulty of tackling a non-convex problem, robust recovery of the depth can be obtained without integrating the normals. By considering the lightings as unknown vector fields, situations usually considered as "hard", like nearby or extended sources, are easily handled. In future work, we will explore joint regularisation of the lightings, the depth and the albedo and study with more care the well-posedness of the inverse problem, the robustness with respect to the numerous parameters, and more efficient numerics.

Fig. 1 .

 1 Fig.1. Robustness to noise and outliers. First and second columns: 4 of the m = 20 images, with σ = 20% and p = 10%; Third column: ground-truth depth map and albedo; Fourth column: results from[START_REF] Ikehata | Robust photometric stereo using sparse regression[END_REF], using a least squares integrator with homogeneous Dirichlet BC; Fifth column: results after L 1 -TV refinement. State-ofthe-art methods provide a rather noisy albedo, and the quality of the final depth map is widely dependent from the choice of a specific integrator. Rather than choosing a robust integrator, as advised e.g. in[START_REF] Durou | Integrating the normal field of a surface in the presence of discontinuities[END_REF], we believe that it makes more sense to use the proposed scheme to directly refine the depth map and the albedo from the images, so as to avoid the propagation of any bias due to normal estimation. Indeed, rather than seeking a piecewise-smooth shape explaining a possibly biased normal field, it makes more sense to look for a piecewise-smooth shape explaining the images themselves.

  Fig.2. Quantitative evaluation of the diffusion process. Left: RMSE in u as a function of the number k of iterations, for several levels of noise (σ, p), with the same reference values of the parameters as in the first test, denoted (α0, θ d,0 , θu,0) in the case of (α, θ d , θu). Right: influence of the model parameter α and of the optimisation parameters θ d and θu on the diffusion process. Convergence is stable but slow if α or θ d is too high (oversmoothed depth map), or if θu is too high (small descent steps). On the contrary, the descent is too fast to be stable if we set a too low value for α or θ d (no regularisation, inducing numerical difficulties), or if θu is too low (large descent steps).

Fig. 3 .

 3 Fig.3. Bias correction by using spatially-varying lightings. Left, from top to bottom: two of the m = 15 input images; associated calibrated lighting intensities; estimated log depth maps using least-squares photometric stereo with a directional model and with the calibrated spatially-varying lightings. Right: relighted views of the estimated sets of 3D-points, using directional (top) or spatially-varying lightings (bottom). Neglecting the spatial variations of the lightings creates a global drift in the reconstruction, which is drastically reduced by considering spatially-varying lightings. This 3D-reconstruction could then be further improved using the proposed L 1 -TV scheme.

Fig. 4 .

 4 Fig. 4. L 1 -TV optimisation with rough initialisation. We show the diffusion process of u (top), ρ (middle) and s 1 (bottom), at k = 0, k = 10, k = 20 and k = 40. The sharp features are first captured by the L 1 data term, and then the TV proximal operators smooth out the residual noise. Realistic estimation can be provided, even with a very rough initialisation.

  Point cloud x, albedo ρ and updated light fields s i , i ∈ [1, m]. Initialisation: Compute the normals and initialise the albedo ρ 0 by per-pixel classical photometric stereo; integrate the normals to initialise u 0

m]; camera parameters; initial light fields s i,0 , i ∈ [1, m]; shape prior u0; model parameters α, β, γ, ǫ; proximal parameters θ d , θρ, θu and θ S ; maximum number of iterations N . Output:

  2. Quantitative evaluation of the diffusion process. Left: RMSE in u as a function of the number k of iterations, for several levels of noise (σ, p), with the same reference values of the parameters as in the first test, denoted (α0, θ d,0 , θu,0) in the case of (α, θ d , θu). Right: influence of the model parameter α and of the optimisation parameters θ d and θu on the diffusion process. Convergence is stable but slow if α or θ d is too high (oversmoothed depth map), or if θu is too high (small descent steps). On the contrary, the descent is too fast to be stable if we set a too low value for α or θ d (no regularisation, inducing numerical difficulties), or if θu is too low (large descent steps).

Acknowledgments. This work is part of a technological transfer between IRIT and the Pixience company, funded by Toulouse Tech Transfer.


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

Jean-Denis Durou

Edge-Preserving Integration of a Normal Field: Weighted Least-squares, TV and L 1 Approaches

Keywords: Integration, Shape-from-gradient, Photometric stereo

We introduce several new functionals, inspired from variational image denoising models, for recovering a piecewise-smooth surface from a dense estimation of its normal field 1 . In the weighted least-squares approach, the non-differentiable elements of the surface are a priori detected so as to weight the least-squares model. To avoid this detection step, we introduce reweighted least-squares for minimising an isotropic TV-like functional, and split-Bregman iterations for L 1 minimisation.

Introduction

Problem Statement. The normal field n of a surface can be estimated by 3Dreconstruction techniques such as photometric stereo [START_REF] Woodham | Photometric method for determining surface orientation from multiple images[END_REF]. To obtain a set of 3D points located on the surface, the estimated normal field must then be integrated into a depth map z, over a subset Ω of the image domain. This second step is crucial in the 3D-reconstruction process, since the accuracy of the recovered surface widely depends on the robustness of integration to noise and outliers.

Let us first recall the equations describing this integration problem, which are similar under both orthographic and perspective projections. In the orthographic case, z is related to the normal field n, for every (x, y) ∈ Ω, through [START_REF] Durou | Integration of a Normal Field without Boundary Condition[END_REF]:

n(x, y) = 1 ∇z(x, y) 2 2 + 1 -∇z(x, y) 1 (1)
where ∇z = [∂ x z, ∂ y z] ⊤ is the gradient of z. Denoting:

p O (x, y) = - n 1 (x, y) n 3 (x, y) , q O (x, y) = - n 2 (x, y) n 3 (x, y) , g O (x, y) = [p O (x, y), q O (x, y)] ⊤ (2)
where n i , i ∈ [START_REF] Agrawal | What is the range of surface reconstructions from a gradient field?[END_REF][START_REF] Bredies | Total generalized variation[END_REF], is the i-th component of n, we obtain from (1) and (2):

∇z(x, y) = g O (x, y) (3)
In the case of perspective projection, we need to know the focal length f of the camera, and to set the origin of image coordinates to the principal point. Introducing the change of variable z = log z, we obtain [START_REF] Durou | Integration of a Normal Field without Boundary Condition[END_REF]:

n(x, y) = 1 ∇z(x, y) 2 2 + 1+∇z(x, y)• 1 f [x, y] ⊤ 2 -∇z(x, y) 1+∇z(x, y)• 1 f [x, y] ⊤ (4)
By setting d(x, y) = xn 1 (x, y) + yn 2 (x, y) + f n 3 (x, y), and:

p P (x, y) = -n 1 (x, y) d(x, y) , q P (x, y) = -n 2 (x, y) d(x, y) , g P (x, y) = [p P (x, y), q P (x, y)] ⊤ [START_REF] Durou | Integrating the normal field of a surface in the presence of discontinuities[END_REF] we get from ( 4) and [START_REF] Durou | Integrating the normal field of a surface in the presence of discontinuities[END_REF], after some algebra:

∇z(x, y) = g P (x, y)

Thus, for both these projection models, one has to solve, in every (x, y) ∈ Ω, the same equation: ∇u(x, y) = g(x, y)

where (u, g) = (z, g O ) in the orthographic case, and (u, g) = (z, g P ) in the perspective one.

Integrating the normal field refers to the process of recovering the unknown u, which will be abusively referred to as "depth map" in the following, from an estimation g = [p, q] ⊤ of its gradient field over Ω. This problem, which has a long history since it dates back to the Dirichlet problem, has given rise to numerous studies in the area of mathematics for imaging, using many different approaches such as Fourier analysis [START_REF] Frankot | A Method for enforcing integrability in shape from shading algorithms[END_REF][START_REF] Simchony | Direct analytical methods for solving Poisson equations in computer vision problems[END_REF], fast marching [START_REF] Galliani | Fast and robust surface normal integration by a discrete eikonal equation[END_REF] or Sylvester equations [START_REF] Harker | Least squares surface reconstruction from measured gradient fields[END_REF][START_REF] Harker | Regularized reconstruction of a surface from its measured gradient field[END_REF]. In this paper, as in many recent works [START_REF] Agrawal | What is the range of surface reconstructions from a gradient field?[END_REF][START_REF] Badri | Robust surface reconstruction via triple sparsity[END_REF][START_REF] Du | Robust surface reconstruction from gradient field using the L1 norm[END_REF][START_REF] Durou | Integrating the normal field of a surface in the presence of discontinuities[END_REF][START_REF] Reddy | Enforcing integrability by error correction using L1-minimization[END_REF], we choose the energy minimisation way, which offers a natural framework for controlling the influence of noise and outliers.

Summary of our Contributions. We focus on the case where solving Eq. ( 7) makes sense only almost everywhere, which happens as soon as the surface to be reconstructed contains edges and depth discontinuities: the gradient ∇u of u cannot be defined on the neighborhood of such non-differentiable elements. In this case, classical least-squares solvers fail (Figure 1) and more robust estimation must be considered. Completing the study proposed in [START_REF] Durou | Integrating the normal field of a surface in the presence of discontinuities[END_REF], we introduce three new functionals inspired by image denoising models, whose minimisation is shown to provide piecewise-smooth surfaces on an arbitrary connected domain Ω. They are based, respectively, on weighted least-squares (WLS), isotropic total variation (TV), and L 1 optimisation.

The rest of this paper is organized as follows. After reviewing in Section 2 the main energy minimisation methods for surface reconstruction from a gradient field, we detail in Section 3 the proposed edge-preserving approaches, which are eventually evaluated on both synthetic and real-world datasets in Section 4.

Fig. 1. Least-squares normal integration. First row, from left to right: ground truth C ∞ depth map u, analytical derivatives p and q corrupted by an additive zero-mean Gaussian noise with standard deviation σ = 5% of g ∞, and least-squares reconstruction [START_REF] Harker | Least squares surface reconstruction from measured gradient fields[END_REF]. Second row: same, for piecewise-C ∞ surface. Noise in the data is successfully handled by least-squares (first row), but discontinuities are smoothed (second row).

2 Related Work

Integrability of a Gradient Field

In the ideal case, g = [p, q] ⊤ is the true gradient of a C 2 function u holding: ∂ yx u = ∂ xy u (Schwarz' theorem). The distance from a gradient field g to an ideal (integrable) field holding ∂ y p = ∂ x q can thus be measured by the integrability term [START_REF] Frankot | A Method for enforcing integrability in shape from shading algorithms[END_REF]:

I(x, y) = |∂ y p(x, y) -∂ x q(x, y)| (8)
which is never null in real-world scenarios, because of noise and of depth discontinuities. In such cases, it makes sense to estimate an approximate solution u of Eq. ( 7) whose gradient ∇u is integrable, rather than to solve Eq. ( 7) exactly. This can be performed efficiently through energy minimisation, by seeking u as the solution of an optimisation problem, lying in an appropriate function space. For instance, if u is sought in L 2 (Ω), integrability of its gradient is implicitely granted (in the presence of discontinuities, the space of functions with bounded variations BV (Ω) should be preferred, so as to allow piecewise-smooth functions). We provide hereafter a brief overview of the main normal integration methods relying on energy minimisation.

Continuous Least-squares Formulation

The most natural energy minimisation approach to solve [START_REF] Frankot | A Method for enforcing integrability in shape from shading algorithms[END_REF] consists in estimating u in a least-squares sense [START_REF] Frankot | A Method for enforcing integrability in shape from shading algorithms[END_REF][START_REF] Simchony | Direct analytical methods for solving Poisson equations in computer vision problems[END_REF][START_REF] Durou | Integration of a Normal Field without Boundary Condition[END_REF][START_REF] Harker | Least squares surface reconstruction from measured gradient fields[END_REF], by introducing the functional:

F LS (u) = Ω ∇u(x, y) -g(x, y) 2 2 dx dy (9)
According to the calculus of variations, minimising this functional is equivalent to solving the associated Euler-Lagrange equation on the interior part

• Ω of Ω: ∆u = ∇• g (10)
which is a Poisson equation (∇• is the divergence operator, which is the adjoint of the gradient, and ∆ = ∇•∇ is the Laplacian operator), along with the natural boundary condition (BC), which is of the Neumann type:

(∇u -g) • µ = 0 (11)
on the boundary ∂Ω of Ω, µ being normal to ∂Ω. Discretising Eqs ( 10) and ( 11) provides a linear system of equations which can be solved in linear time through Fast Fourier Transform (FFT), if Ω is rectangular. Indeed, replacing the natural BC ( 11) by a periodic one, Frankot and Chellappa's well-known algorithm [START_REF] Frankot | A Method for enforcing integrability in shape from shading algorithms[END_REF] recovers the Fourier transform of the depth map analytically, and inverse FFT eventually provides a solution u of [START_REF] Harker | Least squares surface reconstruction from measured gradient fields[END_REF]. This algorithm was extended by Simchony et al. in [START_REF] Simchony | Direct analytical methods for solving Poisson equations in computer vision problems[END_REF] to the natural BC, through the use of Discrete Cosine Transform (DCT).

Discretising the Functional, or the Optimality Conditions?

Instead of discretising the optimality conditions ( 10) and ( 11), the functional (9) itself can be discretised. This is the approach followed in [START_REF] Durou | Integration of a Normal Field without Boundary Condition[END_REF], where it is shown that doing so, no explicit BC is needed (yet, the natural BC is implicitely satisfied). After proper discretisation, a new linear system is obtained, which is solved using Jacobi iterations. Alternatively, Harker and O'Leary show in [START_REF] Harker | Least squares surface reconstruction from measured gradient fields[END_REF] that the discrete least-squares functional can be minimised by solving a Sylvester equation, provided Ω is rectangular (this hypothesis is neither required in [START_REF] Durou | Integration of a Normal Field without Boundary Condition[END_REF], nor in the present paper).

Examples of results obtained using the least-squares solver from [START_REF] Harker | Least squares surface reconstruction from measured gradient fields[END_REF] are shown in the last column of Figure 1. These experiments illustrate the robustness of least-squares against additive Gaussian noise, but also the edge smoothing which occurs using quadratic regularisation. As we shall see later, quadratic regularisation can be improved by introducing weights, or by replacing the squared L 2 norm by a non-differentiable regularisation.

Since the functional ( 9) is convex, discretising the functional or the associated optimality condition should be strictly equivalent, provided the natural BC is enforced. Yet, as noted by Harker and O'Leary in [START_REF] Harker | Regularized reconstruction of a surface from its measured gradient field[END_REF], Poisson-based integration relying on DCT suffers from a bias, due to inconsistent numerical approximations of the gradient ∇u in the discretisation of the natural BC [START_REF] Harker | Regularized reconstruction of a surface from its measured gradient field[END_REF]. The choice of such inconsistent derivatives, as well as a rectangular domain Ω, are actually necessary for obtaining a matrix of the block-Toeplitz type, and thus allowing fast recovery by DCT. Choosing consistent numerical derivatives, or a non-rectangular domain Ω, the structure of this matrix is lost, and the system resulting from the discretisation of the continuous optimality condition must be solved using standard sparse solvers.

In this paper, as in [START_REF] Durou | Integration of a Normal Field without Boundary Condition[END_REF][START_REF] Harker | Least squares surface reconstruction from measured gradient fields[END_REF], we choose to consider discrete functionals so as to avoid dealing with boundary conditions. Rather than relying on special matrix structures [START_REF] Harker | Least squares surface reconstruction from measured gradient fields[END_REF], we use standard solvers for the numerics, allowing us to deal with non-rectangular domains, as in [START_REF] Durou | Integration of a Normal Field without Boundary Condition[END_REF].

Non-quadratic Regularisations

In [START_REF] Harker | Regularized reconstruction of a surface from its measured gradient field[END_REF], Harker and O'Leary extend the method from [START_REF] Harker | Least squares surface reconstruction from measured gradient fields[END_REF] to the case of spectral and Sobolev regularisations, improving the robustness of their method to Gaussian noise. Yet, such regularisations are not adapted to depth discontinuities, since they remain quadratic. In [START_REF] Agrawal | What is the range of surface reconstructions from a gradient field?[END_REF], Agrawal et al. study several functionals having the following general form:

F Ψ (u) = Ω Ψ ( ∇u(x, y) -g(x, y) 2 ) dx dy ( 12
)
where Ψ is chosen so as to reduce the influence of outliers. A numerical study of the discrete versions of several such functionals is presented in [START_REF] Durou | Integrating the normal field of a surface in the presence of discontinuities[END_REF], where the Jacobi iterations used in [START_REF] Durou | Integration of a Normal Field without Boundary Condition[END_REF] are extended to the minimisation of non-convex functionals through semi-implicit schemes. The use of sparse regularisations derived from the L 1 norm has also become an important research direction [START_REF] Du | Robust surface reconstruction from gradient field using the L1 norm[END_REF][START_REF] Reddy | Enforcing integrability by error correction using L1-minimization[END_REF]: we will show how to accelerate such schemes using split-Bregman iterations. Extension to L p minimisation, p < 1, is also presented in [START_REF] Badri | Robust surface reconstruction via triple sparsity[END_REF]. The results are indeed impressive in the presence of very noisy data, but involve setting numerous parameters, which is hardly tractable in real-world applications. Furthermore, since photometric stereo is a technique which is mostly performed inside laboratories, the presence of a huge amount of Gaussian noise in the measurements is very unlikely, and thus greater care is given in this paper to outliers such as discontinuities, which cannot be avoided since they describe the surface itself and not the acquisition procedure. We introduce in the next section several new functionals related to this issue.

3 New functionals

Quadratic Prior

The functional [START_REF] Kimmel | An algebraic multigrid approach for image analysis[END_REF] is not coercive, because of the ambiguity u → u + k, k constant, in the initial equation [START_REF] Frankot | A Method for enforcing integrability in shape from shading algorithms[END_REF]. In the literature, this ambiguity is usually solved a posteriori, for instance by manually setting the mean value of u. In this work, we proceed this way to first compute an approximate solution u 0 through DCT [START_REF] Simchony | Direct analytical methods for solving Poisson equations in computer vision problems[END_REF] (if Ω is not rectangular, g has to be completed with null values, which obviously creates a bias), before introducing u 0 as a quadratic prior to force coercivity.

This prior being biased in the presence of discontinuities and non-rectangular domains, it can be seen as an initial depth map that we want to denoise using an edge-preserving regularisation Φ which shall ensure diffusion along g:

F Φ (u) = Ω Φ (∇u(x, y) -g(x, y)) + λ 2 (u(x, y) -u 0 (x, y)) 2 dx dy ( 13
)
with λ > 0 chosen according to the quality of the approximate solution u 0 .

In this paper, we consider three types of regularisation: Φ = w . 2 2 (weighted least-squares), Φ = . 2 (isotropic TV-like), and Φ = . 1 (L 1 ).

Weighted Least-squares Functional

In a first approach, we assume that it is possible to a priori detect outliers through the evaluation of the integrability term [START_REF] Galliani | Fast and robust surface normal integration by a discrete eikonal equation[END_REF]. This a priori detection is used to weight the influence of discontinuities. Setting Φ(.) = w||.|| 2 2 , where w is a weighting function depending only on the integrability term [START_REF] Galliani | Fast and robust surface normal integration by a discrete eikonal equation[END_REF], and not on u, we obtain the weighted least-squares functional:

F W LS (u) = Ω w(x, y) ∇u(x, y)-g(x, y) 2 2 + λ 2 (u(x, y) -u 0 (x, y)) 2 dx dy (14)
Since we know that the integrability ( 8) is an indicator of the presence of discontinuities (though having null integrability does not imply being smooth: think of a piecewise flat shape), it seems natural to choose for w an integrabilitybased weighting function, which should be a decreasing function of I. To choose effectively the weights, let us use the continuous optimality condition associated with F W LS . Assuming w > 0, and remarking that ∇w w = ∇(log w), we obtain:

∆u + ∇(log w) • (∇u -g) -λ (u -u 0 ) = ∇ • g (15)
Because of the presence of the logarithm, we consider:

w(x, y) = exp(-γ I(x, y) 2 ) ( 16
)
where γ ≥ 0 is a parameter for controlling the weights (γ = 0 corresponds to the standard least-squares formulation). We now discretise u uniformly over a grid (which does not need to be rectangular, unlike in [START_REF] Frankot | A Method for enforcing integrability in shape from shading algorithms[END_REF][START_REF] Simchony | Direct analytical methods for solving Poisson equations in computer vision problems[END_REF][START_REF] Harker | Least squares surface reconstruction from measured gradient fields[END_REF][START_REF] Harker | Regularized reconstruction of a surface from its measured gradient field[END_REF]), with spacing 1, also denoted Ω for convenience. Extending the rationale in [START_REF] Durou | Integration of a Normal Field without Boundary Condition[END_REF] for least-squares functionals, a consistent secondorder accurate discretisation of ( 14) is obtained by first-order forward differences in u, and computation of the forward means of the components p and q of g:

F W LS (u) = (i,j)∈Ω x+ w i,j ∂ + x u i,j -pi,j 2 + (i,j)∈Ω y+ w i,j ∂ + y u i,j -qi,j 2 + λ 2 (i,j)∈Ω u i,j -u i,j 0 2 (17)
where we denote u i,j the value of u at discrete point (i, j), pi,j = p i+1,j +p i,j 2 , qi,j = q i,j+1 +q i,j 2 , ∂ + x u i,j = u i+1,ju i,j , ∂ + y u i,j = u i,j+1u i,j , Ω x+ = {(i, j) ∈ Ω s.t. (i + 1, j) ∈ Ω} and Ω y+ = {(i, j) ∈ Ω s.t. (i, j + 1) ∈ Ω}. The optimality condition in u i,j ∈ Ω reads: χ i+1,j w i,j u i+1,j -u i,j +χ i,j+1 w i,j u i,j+1 -u i,j +χ i-1,j w i-1,j u i-1,j -u i,j +χ i,j-1 w i,j-1 u i,j-1 -u i,jλ 2 u i,j = χ i+1,j w i,j pi,j + χ i,j+1 w i,j qi,j (18)

-χ i-1,j w i-1,j pi-1,j -χ i,j-1 w i,j-1 qi,j-1 - λ 2 u i,j 0
where χ is the characteristic function of Ω. If w is constant and λ = 0, it is easily verified that (18) is a discrete approximation of both the Poisson equation [START_REF] Harker | Least squares surface reconstruction from measured gradient fields[END_REF] and the natural BC [START_REF] Harker | Regularized reconstruction of a surface from its measured gradient field[END_REF].

Stacking the u i,j column-wise in a vector u of size n × 1, where n is the cardinal of Ω, the optimality condition (18) reads as a linear system Au = b, where A is a block-pentadiagonal n × n full-rank matrix with strictly dominant diagonal. We experimentally found that, for relatively small grids (up to 512 × 512), direct sparse solvers provide a fast solution to this system: since A has a small bandwidth (equal to the number of rows in Ω), computation of the sparse product A ⊤ A is very fast, and the normal equation A ⊤ Au = A ⊤ b can be solved through sparse Cholesky factorisation, though it artificially increases the order of points involved in the finite differences, leading to a small additional smoothing (see the numerical results on the peaks dataset in Section 4). Studying more efficient solvers for this problem, as for instance Krylov subspace methods applied to the initial Au = b problem, will be the subject of a future research.

Isotropic TV Functional

The previous approach relies on a priori detection of the discontinuities, so that the corresponding points are "manually" discarded from the equality ( 7). Yet, a priori setting the weights might sometimes be tedious. The weights can also be automatically chosen as a function of ∇u-g 2 [START_REF] Agrawal | What is the range of surface reconstructions from a gradient field?[END_REF][START_REF] Durou | Integrating the normal field of a surface in the presence of discontinuities[END_REF], but the problem cannot be solved directly anymore, and requires an iterative minimisation. We show how to use this idea to minimise a functional resembling the L 2 -TV model [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF].

It is well known in the image processing community that the isotropic total variation (TV) measure T V (u) = Ω ∇u(x, y) 2 dxdy has interesting edgepreserving properties, and tends to favor piecewise-smooth solutions. Considering the discontinuities as the equivalent of edges in image denoising, one would expect the residual ∇ug to be piecewise-smooth as well, with jumps located in discontinuities. This remark invites us to adapt the ROF model [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] to our problem: choosing Φ(.) = . 2 , we obtain from (13) the following functional:

F T V (u) = Ω ∇u(x, y) -g(x, y) 2 + λ 2 (u(x, y) -u 0 (x, y)) 2 dx dy (19) Remarking that . 2 = . 2 . 2 2
. 2 2 , this functional can be minimised through iteratively reweighted least-squares:

       w k (x, y) = ∇u k (x, y) -g(x, y) 2 ∇u k (x, y) -g(x, y) + θ 2 2 , ∀(x, y) ∈ Ω u k+1 = argmin u Ω w k (x, y) ∇u(x, y)-g(x, y) 2 2 + λ 2 (u(x, y)-u 0 (x, y)) 2 dxdy (20)
with u 0 = u 0 , w 0 = 1, and θ > 0, small. The update in u, using Cholesky factorisation, has already been described in Section 3.2. Proceeding so, the normal equations are solved at each iteration, as in [START_REF] Agrawal | What is the range of surface reconstructions from a gradient field?[END_REF]. As a consequence, few iterations are needed, though this might become memory-hungry for large grids. Iterative Jacobi approximations, in the manner of what is proposed in [START_REF] Durou | Integrating the normal field of a surface in the presence of discontinuities[END_REF], would probably offer a less memory-hungry solution. Alternatively, split-Bregman iterations can be considered: we show in the following paragraph how to use such iterations for minimising the anisotropic TV (L 1 ) model.

L 1 Functional

The discontinuities being sparsely distributed in essence, it seems natural to rely on the sparsity enhancing properties of the L 1 norm [START_REF] Du | Robust surface reconstruction from gradient field using the L1 norm[END_REF][START_REF] Reddy | Enforcing integrability by error correction using L1-minimization[END_REF]. Considering the choice Φ(.) = . 1 , we get from (13):

F L 1 (u) = Ω ∇u(x, y) -g(x, y) 1 + λ 2 (u(x, y) -u 0 (x, y)) 2 dx dy (21)
This new functional is still convex, but cannot be minimised through differentiable optimisation. Split-Bregman iterations [START_REF] Goldstein | The Split Bregman method for L1-regularized problems[END_REF] can be considered:

u k+1 = argmin u α 2 d k -(∇u -g) -b k 2 2 + λ 2 u -u 0 2 2 (22)
d k+1 = argmin d d 1 + α 2 d -(∇u k+1 -g) -b k 2 2 (23) b k+1 = b k + (∇u k+1 -g) -d k+1 (24)
where

(d k , b k ) = [d k 1 , d k 2 ] ⊤ , [b k 1 , b k 2 ]
⊤ are auxiliary variables related to the Bregman distance at iteration k. We solve the discrete version of (22) using the same kind of discretisation as in Section 3.2. Yet, unlike in Section 3.2 and 3.3, it is preferable not to solve the problem exactly [START_REF] Goldstein | The Split Bregman method for L1-regularized problems[END_REF], so as to improve convergence properties of the split-Bregman iterations: as advised in the literature, we perform only a few (typically 5) Gauss-Seidel updates at each iteration k. Regarding the basis pursuit problem (23), solution is obtained by shrinkage:

       d k+1 1 = ∂ x u k+1 -p + b k 1 |∂ x u k+1 -p + b k 1 | max ∂ x u k+1 -p + b k 1 , 0 d k+1 2 = ∂ y u k+1 -q + b k 2 |∂ y u k+1 -q + b k 2 | max ∂ y u k+1 -q + b k 2 , 0 (25)
where |.| is the absolute value. We now experimentally compare the proposed schemes on synthetic data, and show results of the L 1 approach on real-world datasets.

Results

Synthetic Data

We first evaluate the performances of the proposed algorithms on synthetic datasets (Figure 2 and Table 1). In each test, a small Gaussian noise with zeromean and standard deviation σ = 0.5% of g ∞ was added to the gradient field, before it was integrated using, respectively, least-squares [START_REF] Simchony | Direct analytical methods for solving Poisson equations in computer vision problems[END_REF][START_REF] Harker | Least squares surface reconstruction from measured gradient fields[END_REF], spectral regularisation [START_REF] Harker | Regularized reconstruction of a surface from its measured gradient field[END_REF], weighted least-squares (λ = 10 -5 , γ = 10), isotropic TV (λ = 10 -5 , θ = 10 -3 ) and L 1 (λ = 10 -4 , α = 0.1). The convergence criterion for the iterative methods was set to a 5.10 -4 relative residual between u k and u k+1 . For fair comparison, the integration constant was changed a posteriori so as to minimise the RMSE between the estimated depth map and the ground truth. The performances of each algorithm are evaluated by Matlab codes running on a I7 laptop at 2.9 Ghz. Fig. 2. Results on synthetic data. We show the ground truth surface (first row), the results using spectral regularisation [START_REF] Harker | Regularized reconstruction of a surface from its measured gradient field[END_REF] (second row), and those using the weighted least-squares (third row), TV (fourth row) and L 1 (fifth row) functionals. L 1 minimisation qualitatively offers the sharper edges, though a staircase effect appears. Weighted least squares provide accurate results for the "Canadian Tent", because the discontinuities correspond to very high integrability values, but it does not perform as well on the "Synthetic Vase", since a part of the discontinuity has null integrability.

Peaks

Canadian Tent Synthetic Vase (512 × 512) (256 × 256) (320 × 320) Least-squares (DCT) [START_REF] Simchony | Direct analytical methods for solving Poisson equations in computer vision problems[END_REF] 0.30 (0.05 s) 10.76 (0.02 s) 4.55 (0.03 s) Least-squares (Sylvester) [START_REF] Harker | Least squares surface reconstruction from measured gradient fields[END_REF] 0.14 (0.84 s) 10.76 (0.32 s) 4.56 (0.46 s) Spectral regularisation [START_REF] Harker | Regularized reconstruction of a surface from its measured gradient field[END_REF] 0. 1. RMSE (in pixels) between the ground truth depth maps and those recovered using three state-of-the-art algorithms and our three new ones. The "Peaks" depth map being C ∞ , all methods succeed at recovering accurate results for this dataset. Since we solve the normal equation (by means of Cholesky factorisation) in the WLS and TV approaches, additional smoothing is introduced, and thus these methods perform a little better than L 1 in this test.

Applications on Real Data

Photometric Stereo. The proposed split-Bregman scheme (Section 3.4) was applied to real-world gradient fields (the other proposed schemes provide comparable results), obtained by applying the photometric stereo technique [START_REF] Woodham | Photometric method for determining surface orientation from multiple images[END_REF] to the "Scholar"2 and to the "Beethoven"3 datasets (Figure 3). To emphasize the discontinuity-preserving properties of the scheme, as well as the staircasing effect appearing on large flat areas, we applied the method on the whole rectangular domain, rather than manually segmenting Ω.

It should be noted that a staircasing effect occurs on the background. This effect is well known and studied in the context of image denoising: adapting the total generalised variations schemes (TGV) [START_REF] Bredies | Total generalized variation[END_REF], we could probably get rid of it.

Yet, staircasing seems to affect only the background, and should thus not be considered as a really damaging effect, since our method is able to deal with non-trivial integration domains (which is not the case in many algorithms [START_REF] Frankot | A Method for enforcing integrability in shape from shading algorithms[END_REF][START_REF] Simchony | Direct analytical methods for solving Poisson equations in computer vision problems[END_REF][START_REF] Harker | Least squares surface reconstruction from measured gradient fields[END_REF][START_REF] Harker | Regularized reconstruction of a surface from its measured gradient field[END_REF]): staircase-free 3D-reconstructions can be obtained by manually segmenting the reconstruction domain (Figure 4). We studied weighted least-squares, TV and L 1 functionals in the context of normal field integration, and provided efficient numerics for minimising these functionals, through sparse Cholesky factorisation, reweighted least-squares and split-Bregman iterations, respectively. We experimentally showed that these functionals provide sharp depth maps, and demonstrated how to use them in the context of photometric stereo and surface edition. In future work, we plan to accelerate the iterative schemes by multigrid techniques [START_REF] Kimmel | An algebraic multigrid approach for image analysis[END_REF], so as to allow realtime surface reconstruction and edition, and to more deeply study the staircase effect appearing with L 1 minimisation. We believe that introducing higher order regularisation terms, in the spirit of the total generalised variation regularisation (TGV) [START_REF] Bredies | Total generalized variation[END_REF], will annihilate this effect, while providing a higher order of accuracy.

Fig. 3 .

 3 Fig. 3. Photometric stereo. A scene is captured from the same point of view, but under different lightings (top row), so that the normal field can be revealed using photometric stereo. By integrating this normal field (split-Bregman iterations), we obtain the surfaces on the bottom row. The staircase effect is clearly visible in these examples, though it only affects the background.

Sample codes for testing the proposed methods can be found on http://ubee. enseeiht.fr/photometricstereo/

http://vision.seas.harvard.edu/qsfs/Data.html

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

Surface Edition. In order to further illustrate the iterative normal integration through split-Bregman iterations, we consider a surface edition problem, consisting in inserting m small objects, whose gradient fields g i , i = 1 . . . m, are known over Ω i ⊂ Ω, into a larger object represented by its gradient field g 0 over Ω and the corresponding least-squares depth map u 0 , while preserving thin details. Setting Ω 0 = Ω \ ∪ m i=1 Ω i , this reads as the minimisation of:

(26) which is an extension of the Poisson image editing problem [START_REF] Pérez | Poisson image editing[END_REF]. Functionals (21) and ( 26) are the same, provided that g = m i=0 χ Ω i g i , where χ Ω i is the characteristic function of Ω i . Now, we merge both the gradient fields g 0 of the "Scholar" dataset (the reconstruction domain Ω is set to the non-rectangular domain of this dataset), and g 1 of the "Beethoven" dataset, so as to replace a small area Ω 1 ⊂ Ω of the reconstructed "Scholar" surface, by Beethoven's bust. In addition, we would like to remove some details inside another domain Ω 2 ⊂ Ω (Figure 4). To this purpose, we choose g 2 = 0: this will perform TV-"inpainting" inside Ω 2 . Denoting Ω 0 = Ω \ Ω 1 ∪ Ω 2 , we form the gradient field g = 2 i=0 χ Ω i g i , and apply the proposed split-Bregman scheme. As shown in Figure 4, a detail-preserving blending of the statues is obtained, while removing the details inside Ω 2 . Fig. 4. Surface edition. In the top-left figure, the colored part is Ω, the red part is Ω 0 ("Scholar" gradient field), the yellow one is Ω 1 ("Beethoven" gradient field), and the purple one is Ω 2 (inpainting area). We show the diffusion process at iterations 0, 10 (fine details begin to appear in the "hair"), 50 (the details inside the inpainted area disappear), 200 and 1000 (stable). Apart from the initialisation u0 using the DCT solver [START_REF] Simchony | Direct analytical methods for solving Poisson equations in computer vision problems[END_REF], which required to add the background to Ω with null values of the gradient, the background was removed from the reconstruction domain: this considerably improves the boundaries of the surface, proving the importance of considering integration schemes able to deal with non trivial domains.


Franck Ravat
email: ravat@irit.fr

Jiefu Song
email: song@irit.fr

Olivier Teste
email: teste@irit.fr

Analysis Framework for Reduced Data Warehouse

Our aim is to define a framework supporting analysis in MDW with reductions. Firstly, we describe a modeling solution for reduced MDW. A schema of reduced MDW is composed of states. Each state is defined as a star schema composed of one fact and its related dimensions valid for a certain period of time. Secondly, we present a multi-state analysis framework. Extensions of classical drilldown and rollup operators are defined to support multi-states analyses. Finally we present a prototype of our framework aiming to prove the feasibility of concept. By implementing our extended operators, the prototype automatically generates appropriate SQL queries over metadata and reduced data.

Introduction

Nowadays, Multidimensional Data Warehouse (MDW) is a widely used component in decision support systems. A MDW schema is based on facts (analysis subjects) and dimensions (analysis axes). The facts contain analysis indicators while the dimensions organize analysis parameters according to their level on hierarchies: from the minimal (most detailed) granularity to the maximal (most general) granularity.

In a MDW, data are stored permanently and new data is steadily added. As result, a MDW stores a huge volume of data in which the analyst may be lost during her/his analyses. On the other hand, the relevance of MDW data decreases with age: detailed information is generally considered essential for recent data [START_REF] Skyt | Specification-based data reduction in dimensional data warehouses[END_REF], while more aggregated information can usually satisfy the need of analysis over older data. For instance, an analyst may have interest in analyzing published news by subthemes for the last four years. However, as most of today's subthemes did not exist before, the subtheme granularity level may be proved useless for an older period. As a result the analyst may have no more interest in analyzing published news by subtheme over the last ten years but by a higher and more stable granularity level, such as news' theme.

Facing large volumes of data among which a great amount of inadequate data are found, our aim is to both increase the efficiency of analysis and facilitate the analysts' task. To this end, we provide a conceptual MDW model which provides only pertinent data over time. Meanwhile we also develop a framework compatible with our conceptual MDW model in order to provide powerful analysis tools to decision-makers.

-The MDW model with reduction increases analysis efficiency by allowing analysts to remove useless temporal granularity levels according to their needs. As detailed information loses its value over time, we intend to implement selective deletion at low granularity levels. -The compatible framework permits to manage reduced data and model analysis process by including analysis operators and a graphical interface.

This paper is composed of the following sections. Section 2 describes a state of the art of data reduction. Section 3 presents preliminary concepts of reduced MDW with the help of a case study. Section 4 describes our analysis framework compatible with reduced MDW by emphasizing modeling principles of metamodel. Section 5 focuses on extended analysis operators. Section 6 presents a prototype showing the feasibility of concepts.

Related work

Reducing data allows us to both decrease the quantity of irrelevant data in decision making and increase future analysis quality [START_REF] Udo | Hybrid Data Reduction Technique for Classification of Transaction Data[END_REF] . In the context of decision support, data reduction is a technique originally used in the field of data mining [START_REF] Okun | Unsupervised data reduction[END_REF] [START_REF] Udo | Hybrid Data Reduction Technique for Classification of Transaction Data[END_REF].

In the data warehouse context, [START_REF] Garcia-Molina | Expiring Data in a Warehouse[END_REF] were the first to propose solutions for data deletion. They study data expiration in materialized views so that they are not affected but maintained after updates with the help of a set of standard predefined views. No discussion about carrying out analyses in data warehouse after data expiration can be found in this work.

In the multidimensional approaches, [START_REF] Chen | Olaping stream data: Is it feasible?[END_REF] propose an architecture allowing the integration of data streams into a MDW by reduce their size. The size reducing process is predefined and automatically executed by partially aggregating the data cube; it makes sure the detailed information is only available for a certain period of time. But this work only focuses on the fact table, and no analysis support component is included in proposed architecture. [START_REF] Skyt | Specification-based data reduction in dimensional data warehouses[END_REF] present a technique for progressive data aggregation of a fact. This study intends to specify data aggregation criteria of a fact due to higher levels of dimensions. Although the authors propose techniques to query reduced multidimensional data, they fail to provide metamodel permitting to manage data reduction in MDW. [START_REF] Kimball | The data warehouse toolkit: the complete guide to dimensional modeling[END_REF] define the concept of Slowly Changing Dimensions (SCDs) and indicate that data may change within a dimension even though it occurs less frequently than in a fact. They propose three basic modeling solutions for managing dimensional data changes, namely overwrite old data, create new record for each change and keep data changes as alternative values in MDW. [START_REF] Golfarelli | A survey on temporal data warehousing[END_REF] point out that not only data but also MDW schema can change over time according to user's requirements. SCD, however, does not provide solution for handling schema changes and it does not take user's need into account.

In [START_REF] Iftikhar | A rule-based tool for gradual granular data aggregation[END_REF], a gradual data aggregation solution based on conception, implementation and evaluation is proposed. This solution is based on a table containing different temporal granularities: second, minute, hour, month and year. Unfortunately this work does not discuss possibilities of carrying out analysis over gradual aggregated data.

The previous works only focuses on reduction of fact table or data changes within dimensions. [START_REF] Iftikhar | Using a Time Granularity Table for Gradual lar Data Aggregation[END_REF] and [START_REF] Iftikhar | A rule-based tool for gradual granular data aggregation[END_REF] -We propose a solution for handling the complete MDW schema changes by generalizing the mechanism of reduction. In consequence all of the dimensions as well as the fact are susceptible to sustain reductions to different granularity levels so as to fully satisfy analyst's needs. The information judged useless is aggregated and then deleted from MDW in order to provide only necessary data for analysis. -We face issues of carrying out analysis in MDW whose schema changes over time. We provide decision-makers with an analysis-support framework applicable to MDW with reduction. The framework contains allows managing both reduced MDW and analysis process.

3 Preliminary concepts of reduced MDW

We will firstly describe a case study of data reduction in MDW that fulfills decision-makers' needs. This case study intends to give a first glance at selective deletion of data in MDW as well as data model notations. Then we define a formal presentation for reduced MDW. At the end of the section we describe a set of reduction operators permitting to define reduced MDW.

Case study

This case study shows a complete multidimensional schema progression that fulfills the analyst's needs. A MDW populated by the RSS streams allows decision-makers to analyze the number of published news from her/his favorite websites. Containing over ten million tuples in the fact table, this MDW face the first important problem of performance while a decisionmaker carries out analyses. Moreover, most of the old detailed data become obsolete nowadays, they are no more used in analysis process and should be deleted according to user' needs. More precisely, a decision-maker expresses her/his needs as followed: (a) during the last four years, news analysis is carried out with reference to lowest levels of granularity (subtheme, city and publication date); (b) in the previous period from 2000 to 2010, analyses are summarized according to news' theme, country mentioned in the news and month of publication because no daily analysis referring to subtheme and city is required; (c) before 2000, only aggregated information about published news by quarter and by continent makes sense.

The following three figures represent the conceptual multidimensional schemas fulfilling user's needs. Each schema is based on star schemas introduced in [START_REF] Golfarelli | Conceptual design of data warehouses from e/r schemes[END_REF]. A star schema is based on a subject of analysis (fact) related to different dimensions. Each fact is composed of one or more indicators. For instance, in figure 1 the fact named "FNews" contains one indicator: number of published news (NBN). A dimension models an analysis axis; it represents information according to which subjects of analysis are to be dealt with. For instance, the "FNews" fact is associated to 3 dimensions: DTheme, DGeography and DTimes.

Dimension attributes (also called parameters or levels) are organized according to one or more hierarchies such as HTHM on dimension DTheme.

Concepts

In order to manage progressive schema changes, we model a reduced MDW with a set of star schemas called states. The current state corresponds to the present status of the MDW containing the most detailed information, while past states correspond to a succession of reduced states over time in which information is aggregated.

-A MDW is defined by S = (n S ; E; Map), where n S ∈ N is the name of the MDW; E = E 1 ; ...; E n is a set of states composing the MDW; Reduction function Map : E → E|Map(E k ) = E k+1 defines the state named E k+1 obtained by the reduction of E k . -Each state Ei = (F i ; D i ; T i ) is a star schema defined for a temporal period, where F i ∈ F is a fact representing a subject of analysis; D i = {D times ; D 1 ; ...; D m } ⊆ D is a set of dimensions associated to the fact with necessarily a temporal dimension denoted D times .

To define the T i temporal interval, we adopt a linear and discrete time model approaching time in granular way through time observation units [START_REF] Wang | Logical design for temporal databases with multiple granularities[END_REF]. A temporal interval T k of state E k is defined by a couple of instants. These instants can be fixed (temporal grains such as the year of 1990) or dynamic (defined with the instant t now ). -A fact denoted F i is defined by

F i = (n F i , M F i ), where n F i is the fact name; M F i is a set of measures or indicators, ∀i, j|i < j → M j ⊆ M i . A dimension denoted D i is defined by (n Di , A Di , H Di )
, where n Di is the dimension name; A Di is the set of the attributes of the dimension;

H Di is a set of hierarchies, ∀i, j|i < j → D j ⊆ D i ∧ H j ⊆ H i ∧ A j ⊆ A i . -A hierarchy, denoted H j (abusive notation ofH Di j
) is defined by (n Hj , P Hj , ≺ Hj , W eak Hj ), where n Hj is the hierarchy name; P Hj is a set of attributes called parameters; ≺ Hj is an antisymmetric and transitive binary relation between parameters; W eak Hj is an application that associates to each parameter a set of dimension attributes, called weak attributes. Hierarchies organize the attributes of a dimension, from the finest graduation (root parameter denoted ID Di ) to the most general graduation (extremity parameter denoted ALL Di ). Thus, a hierarchy defines the valid navigation paths on an analysis axis.

-Analysis results are presented in forms of multidimensional table, denoted M T i , ∀i ∈ [1.

.n], which is defined by

(S i ; Ax i ; R i ; I i ), where S i ∈ F is the analysis subject; Ax i = {ax i 1 , ax i 2 , [, ax i 3 [, ax i 4 , . . .]]
} is the set of analysis axes currently presented, among which ax i 1 and ax i 2 are the horizontal and vertical displayed axes respectively; R i =< pred 1 , pred 2 , . . . > is a set of selection predicates filtering displayed analysis results; I i ⊆ D T IM ES is a set of temporal interval representing the validation period of M T i .

Multi-states analysis framework

Data reduction brings changes to not only MDW's schema but also its instances. In regard to schema changes, a MDW is no more modeled as single static schema but a set of schemas (states) over time. As for instance changes, a reduced MDW no longer contains all detailed information over time. Judged useless by analysts, certain detailed information in recent states would be aggregated and then deleted. As consequence, dimensions and facts in recent and older states certainly have different instances, as they contain different attributes and measures respectively.

Due to the schema and instances changes, existing metamodel for MDW as well as analysis operators is no more applicable to reduced MDW. In order to manage different states in MDW and provide decision-makers with efficient analysis tools, we propose an analysis framework compatible with reduced MDW. In this section, we first present the architecture of our framework and the roles of its mains components. Then we detail Data Management parts along with a deeper discussion about how reduced MDW can be managed through our framework. Details about Analysis Engine can be found in section 5.

Architecture of multi-states analysis framework

The multi-states analysis framework is composed of three components, namely Data Management, Analysis Engine and Interactive Restitution (cf. figure 4). Each part has specific roles and interacts with others. This multi-states analysis framework guaranties the transparency of data reduction in MDW. Decision-makers carry out analysis via graphical implementation of algebraic operators (arrow tagged "1" in figure 4) and then receive a global analysis result (arrow tagged "11" in figure 4). No knowledge about schema and instance evolution of MDW is required for effectuating analysis via our framework. This is thanks to the Data Management and Analysis Engine parts of our framework which adapt themselves to users' needs.

Data management

The The ternary association in our metamodel permits to organize attributes according to their level on a hierarchy of a dimension, which corresponds to concept hierarchy H j = (n Hj , P Hj , ≺ Hj , W eak Hj ).

Analysis

One of the core components of our framework is Analysis Engine. In this section we present more details about multi-states analysis processing steps. Detailed studies about emblematic and the most used operators (drilldown and rollup) are also effectuated. To explain how elementary analysis operations are carried out, we present execution algorithms for algebraic operators.

Analysis processing

To facilitate the decision-maker's tasks, she/he only interactively selects the MT the displayed measures of a fact, the displayed attributes of dimensions as well as the temporal interval of the analysis. The system converts the selection into corresponding analysis operators. All operators supporting multi-states analysis are processed according to the five steps showed in figure 6. More precisely, through a set of temporal intervals chosen by decisionmakers, the analysis operator firstly determines the states in which analysis is carried out before adjusting input parameters if necessary. Then it splits up into several classical operators, each one of them is applied to a single concerned state along with input parameters eventually adjusted. Next each partial operator is translated and executed independently to get partial results. At last all partial results are combined in order to return a unique global result to decision-makers.

FIG. 6: Analysis processing steps of multi-states analysis operators

Our proposition of analysis processing differs from the one in [START_REF] Morzy | On querying versions of multiversion data warehouse[END_REF]. Firstly, [START_REF] Morzy | On querying versions of multiversion data warehouse[END_REF] involves only issues at querying level, no discussion about analysis operators can be found. Secondly, decision-makers should provide the exact set of schemas in which analysis is carried out, which requires her/him to have a profound knowledge about schema evolution in MDW. At last, analysis results are possible to be presented separately, which makes the interaction with data more complex.

Conversely, our proposition of analysis processing has multiple advantages. First of all, it depends on formally defined analysis operators which are extensions of classical operators. Existing analysis tools and languages are fully reusable to the extended operator. Secondly, data reduction in MDW is entirely transparent to decision-makers. While carrying out analysis, decision-makers do not necessarily need to know schema evolutions in MDW. The multi-states operators systematically search concerned states and adapt themselves if necessary. Moreover, multi-states analysis operators provide a single analysis result even if information may come from several states. No information about states in MDW is presented to decision-makers, which reinforces the transparency of data reduction.

Extended analysis operators

Ravat et al. (2008) define a set of generic algebraic analysis operators based on MDW with one stable schema over time. Even though these analysis sustain flexible and adaptable to different classical modeling solutions (star and constellation models), they become incompatible with reduced MDW. Firstly, classical analysis operators are conceived to manipulate one unique schema rather than a set of schemas evolving over time. Secondly, classical analysis operators do not provide solutions for handling heterogeneities of instances in different states. Conceived to modify analysis precision, the classical drilldown and rollup operators are highly affected by data reduction in MDW. An analysis involving an inexistent parameter, for instance, is theoretically impossible according to classical analysis operators.

Facing to these issues, we propose two extended analysis operators: Drilldown multi-states and Rollup multi-states , in order to support multi-states analyses in reduced MDW. They both take a set of temporal intervals, a multidimensional table currently displayed, an analysis dimension and a parameter as input. A new multidimensional table is produced as output containing information at lower or higher granularity level. As data reduction is completely transparent, decision-makers have no information about which dimension D i in state E i to choose to carry out analyses. For the sake of simplicity, we allow decision-maker to simply specify the name of dimension n Di in which analyses are carried out. The analysis operators adapt themselves by finding the corresponding dimensions D i in state E i while executing. In order to distinguish dimensions D i in state E i (cf. section 3.2) from dimension chosen by decisionmakers, in the following paper we add "_E i " as suffix to dimension in state E i , while the name of dimension chosen by decision-makers is without suffix.

Drilldown

S K+1 = S K -Analysis axesAx K+1 ⊆ D, ∀ax K+1 j ∈ Ax K+1 |ax K+1 j = (D K+1 j , H K+1 j , P K+1 j ) such as -Dimension D K+1 j = D K j , -Hierarchy H K+1 j = H K j , -Parameter P K+1 j ⊆ P ( H K+1 j )|(D K+1 j = D i → P K+1 j =< P inf > ∨ < P λ 1 >) ∧ (D K+1 j = D i → P K+1 j =< P D K j 1 , P D K j 2 , . . . >) -Selection predicates on dimensions and/or fact R K+1 = R K -Validation period I k+1 = I ∨ I k+1 = I k 1. Adjusted parameter P λ ∈ P H K+1 j |P inf ≺ H K j P λ
TAB. 1: Algebraic multi-states drilldown operator TheDrilldown multi-states operator permits to display information at a granularity level on currently displayed dimension in several states. As certain parameters of low granularity exist only in most recent states but not in former ones, it is very likely that some input parameters are not present in all states involved in analysis, especially for Drilldown Find the subset of states

E i → ∀E j ∈ E i |I Ej ∈ I ∨ I M T k -Adaptation 7 Let P Drilldown = P 8 Let r = F ALSE 9 While ALL D ≺ H actual P Drilldown ∧ r = F ALSE 10 If ∀E j ∈ E i |P Drilldown ∈ A D
)into query Q -Translation 22 M T Ej K+1 = Result of query Q -Execution 23 M T K+1 = M T K+1 M T Ej K+1 -Integration 24
End for 25

End if 26 End if For instance, if a decision-maker carries out an analysis to City level in "E1" and "E2" states, the initial Drilldown multi-states ([2000,2014], M T k , DGeography, City) operator will be handled as follows: (a) the first step Adaptation replaces inexistent parameter by the nearest superior parameter. In this case, City is replaced by Country which is the nearest superior common parameter on hierarchy HGEO of dimensions DGeography_E1 and DGeography_E2. We obtain by a new operator with an adjusted parameter Drilldown multi-states ( [2000, 2014], M T k , DGeography, Country); (b) the second step Decomposition decomposes the adjusted operator into several classical drilldown operators. Each one of them is then applied to a cube containing one single state deducted from the set of intervals; (c) next we find the classical Translation phase which transforms each drilldown operator into an independent SQL query; (d) the following phase Execution executes independently each SQL query to obtain a set of partial analysis results in forms of cube; (e) at last Integration phase gathers all individual results obtained from each state in order to form one single analysis result. The Rollup multi-states operator consists in moving from finer granularity data to coarser granularity data on a currently displayed dimension in several states (cf. table 2).

Rollup

Rollup multi-states operator reveals data at a higher granularity level. As currently displayed granularity levels exist in all concerned states, chosen parameter denoted P sup presents no doubt in all concerned states. For example, if the decision-makers wants to roll up analysis level from Country to Continent after the previous analysis, she/he can simply carry out the following operator: Rollup multi-states ( [2000, 2014] Find the subset of states

E i → ∀E j ∈ E i |I Ej ∈ I ∨ I M T k 7 For E j ∈ E i -Decomposition 8 Let M T Ej K be the part of MT in states E j 9 Translate Rollup(M T Ej K , D, P )into query Q -Translation 10 M T Ej K+1 = Result of query Q -Execution 11 M T K+1 = M T K+1 M T Ej K+1 -Integration 12
End for 13 End if

Implementation

Based on our proposed framework, we develop a prototype in Java. The prototype implements the reduced MDW of our case study. As our previous work [START_REF] Atigui | Reducing multidimensional data[END_REF] shows that queries are more efficiently computed within a reduced Star Schema in DBMS Oracle 11g, we chose to implement our case study with the R-OLAP schema presented in figure 7. As a first step, the aim of our prototype is to demonstrate the feasibility of carrying out multi-states analysis in reduced MDW. Issues about execution efficiency of multi-states analysis operators will be included in our future work. Moreover, we choose to implement in the first place the reduced MDW in a R-OLAP environment, because R-OLAP suffers from poorer query execution performance than M-OLAP and H-OLAP [START_REF] Vassiliadis | A survey of logical models for olap databases[END_REF], especially in an era of massively abundant data. Thus, one of the main focuses of our future work is to evaluate analysis efficiency in other OLAP environments, such as M-OLAP and H-OLAP. Now we will illustrate how multi-states analysis operators are transformed into SQL queries over metadata and reduced data. A decision-maker is analyzing published news in the world by month from 2000 to 2014 (states "E1" and "E2" in reduced MDW). Wishing to consult the number of published news by month and by continents to get more detailed information, she/he slides the cursor on line to Continent level. The framework detects it concerns a multi-states drilldown operation from ALL_G to Continent on dimensions DGeography_E1 and DGeog-raphy_E2. After executing a Drilldown multi-states operator, two queries are generated by Analysis Engine of framework.

-By translating algebraic Drilldown multi-states operator, the prototype generates firstly the query applicable to metamodel (cf. left part in figure 8). It interrogates the metamodel in order to find states in which the drilldown operation is carried out and verify if the chosen parameter is available in concerned states. In this case parameter Continent presents in both "E1" and "E2" states. This paper resides within the field MDW. Our first objective is to specify reduced multidimensional schema over time in order to store only the useful data for decision support according to the needs of analysts. The second objective is to provide a framework allowing analysts to carry out multi-states analysis in reduced MDW.

Firstly, we define a MDW model which allows us to specify MDW schemata composed of a set of states varying over time. Each state consists of a star schema valid for a certain period of time. Secondly, we present a generic framework permits to manage both reduced and unreduced MDW through a metamodel and support multi-states analysis with the help of extended analysis operators. By instantiating the metamodel, we explain how reduced MDW can be managed by our framework. As regard to multi-states analysis operators, we present an algebraic form followed by an execution algorithm in order to show how analysis results are produced. One of the main advantages of our analysis framework is the transparency of data reduction: decision-makers can freely carry out analysis without the need for knowing schema evolution in MDW. Finally, we implement our multi-states analysis framework in order to show the feasibility of proposed concepts. We show how multi-states analysis operators are executed by giving concrete examples of automatically generated SQL queries.

In the future, we intend to integrate more analysis operators in a short term, such as Dice, Slice, etc. Moreover, even though our previous work has shown that queries were more efficiently computed in reduced R-OLAP MDW Atigui et al. (2014), the efficiency of multi-states analysis in M-OLAP and H-OLAP environments still remains to be evaluated. Our long term goals are to study influence of data reduction over pre-aggregated data.

FIG. 1 :

 1 FIG. 1: MDW schema valid from 2010 to 2014

FIG. 4 :

 4 FIG. 4: Main components of multi-states analysis frameworks

  FIG. 5: UML class diagram metamodel of the reduced MDW

  multi-states operator Drilldown multi-states ([I] ; M T k , D i ; P inf )= M T K+1 ; Input -I ⊆ D T IM ES : a set of optional temporal intervals on dimension D T IM ES -M T k : multidimensional table currently displayed -D i : analysis axis currently displayed -P inf : chosen parameter on dimension D i Output -M T K+1 = (S K+1 ; Ax K+1 ; R K+1 ; I k+1 ) analysis result MT such as -Analysis subject

∈

  multi-states operator Rollup multi-states ([I] ; M T k , D i ; P inf )= M T K+1 ; Input -I ⊆ D T IM ES : a set of optional temporal intervals on dimension D T IM ES -M T k : multidimensional table currently displayed -D i : analysis axis currently displayed -P sup : chosen parameter on dimension D i Output -M T K+1 = (S K+1 ; Ax K+1 ; R K+1 ; I k+1 ) analysis result MT such as -Analysis subject S K+1 = S K -Analysis axesAx K+1 ⊆ D, ∀ax K+1 j on dimensions and/or fact R K+1 = R K -Validation period I k+1 = I ∨ I k+1 = I k TAB. 2: Algebraic multi-states rollup operator

  FIG. 7: R-OLAP schema of case study

  use a temporal table for gradual data reduction. Analysts' needs are ignored in data reduction process. None of the previous work fully support carrying out analysis in reduced MDW. Facing to these issues, our goal is more ambitious as it aims to meet several challenges mentioned in[START_REF] Golfarelli | A survey on temporal data warehousing[END_REF]:

  multi-states operator. For example, from 2000 to 2014, hierarchy HGeo on dimensions DGeography_E1 and DGeography_E2 can sustain a drilldown operation to Continent and Country levels. But a classical drilldown until to City level is impossible because the parameter City exists no longer in "E2" state after reduction. To handle this issue, we propose to augment parameter's granularity level until finding the first common parameter among all involved states(cf. table 1). The execution algorithm of Drilldownmulti-states is as follows. Drilldown multi-states ([I] ; M T k , D i ; P inf ) Input: Set of temporal intervals I, displayed multidimensional table M T k , displayed dimension D, parameter P. Output: new multidimensional table M T k+1 1 Let H actual be the actually displayed hierarchy 2 Let P actual be the actually displayed parameter 3 If P actual ≺ H actual P ∨ ALL D ≺ H actual P then

	Algorithm 1: 4 Impossible operation
	5 Else
	6

  , M T k+1 , DGeography, Continent). No special treatment is needed in Adaptation step (cf. algorithm 2). Rollup multi-states ([I] ; M T k , D i ; P inf ) Input: Set of temporal I, displayed multidimensional table M T k , displayed dimension D, parameter P. Output: new multidimensional table M T k+1 1 Let H actual be the actually displayed hierarchy 2 Let P actual be the actually displayed parameter 3 If P ≺ H actual P actual then

	Algorithm 2: 4 Impossible operation
	5 Else
	6