Marco Bevilacqua

Jean-François Aujol

Pierre Biasutti

Mathieu Brédif

Aurélie Bugeau

Joint Inpainting of Depth and Reflectance with Visibility Estimation $

Keywords: Inpainting, Total Variation, Depth Maps, Lidar, Reflectance, Point Cloud, Visibility

This paper presents a novel strategy to generate, from 3-D lidar measures, dense depth and reflectance images coherent with given color images. It also estimates for each pixel of the input images a visibility attribute. 3-D lidar measures carry multiple information, e.g. relative distances to the sensor (from which we can compute depths) and reflectances. When projecting a lidar point cloud onto a reference image plane, we generally obtain sparse images, due to undersampling. Moreover, lidar and image sensor positions typically differ during acquisition; therefore points belonging to objects that are hidden from the image view point might appear in the lidar images. The proposed algorithm estimates the complete depth and reflectance images, while concurrently excluding those hidden points. It consists in solving a joint (depth and reflectance) variational image inpainting problem, with an extra variable to concurrently estimate handling the selection of visible points. As regularizers, two coupled total variation terms are included to match, two by two, the depth, reflectance, and color image gradients. We compare our algorithm with other image-guided depth upsampling methods, and show that, when dealing with real data, it produces better inpainted images, by solving the visibility issue.

Introduction

Image-based 3D reconstruction of static and dynamic scenes [START_REF] Herbort | An introduction to image-based 3D surface reconstruction and a survey of photometric stereo methods[END_REF][START_REF] Seitz | A Comparison and Evaluation of Multi-View Stereo Reconstruction Algorithms[END_REF][START_REF] Stoykova | 3-D time-varying scene capture technologies-A survey[END_REF] is one of the main challenges in computer vision nowadays. In the recent years many efforts have been made to elaborate configurations and approaches, possibly requiring the employment of multiple sensors, with the final goal of generating plausible and detailed 3D models of scenes. To this end, typical optical cameras are often combined with non-visual sensors. The intermediate outputs of these hybrid systems, prior to the final scene rendering, are in general depth or depth+color images (RGB-D). Among the non-visual sensors, we can find Time-of-Flight (ToF) cameras [START_REF] Kolb | Time-of-Flight Cameras in Computer Graphics[END_REF], which acquire low-resolution co-registered depth and color images at a cheap cost, and the famous Kinect [START_REF] Zhang | Microsoft Kinect sensor and its effect[END_REF], capable to extract depth information by exploiting structural light. Another possibility is represented by lidar devices, which are used in a variety of applications and provide as output point clouds with measures of distance and reflectivity of the sensed surfaces. This work lies in the context described and is particularly driven by the exploitation of data acquired by Mobile Mapping Systems (MMS), such as [START_REF] Paparoditis | Stereopolis II: A multi-purpose and multi-sensor 3D mobile mapping system for street visualisation and 3D metrology[END_REF]. MMS systems are vehicles equipped with highresolution cameras and at least one lidar sensor: their contained dimensions allow them to be driven through regular streets and acquire data of urban scenes. The data acquired is a set of calibrated and geolocated images, together with coherent lidar point clouds. The interest towards them comes from the possibility of having available, at a relatively small processing cost, the combination of depth and color information, without having to perform explicit (error-prone) reconstructions. Having a good depth estimate at each pixel, for example, would enable the possibility to perform depth-imagebased rendering algorithms, e.g. [START_REF] Chen | Efficient Depth Image Based Rendering with Edge Dependent Depth Filter and Interpolation[END_REF][START_REF] Schmeing | Depth Image Based Rendering[END_REF][START_REF] Zinger | Free-viewpoint depth image based rendering[END_REF]. Similarly, the availability of depth information allows the insertion of virtual elements into the image, such as pedestrians or vehicles generated by a traffic simulation [START_REF] Brédif | Image-Based Rendering of LOD1 3D City Models for traffic-augmented Immersive Street-view Navigation[END_REF]. While MMS data sets do not include directly depth images aligned with the available color images, it is easy, by exploiting the known geometry, to project the lidar point clouds onto each image. This operation produces initial depth images, which present three main issues (see Figure 1, where three parts of an input depth image are shown, together with the corresponding image parts). 1. Undersampling: since lidar and image acquisitions are deeply different in terms of geometry and characteristics, the resulting depth images turn to be irregular. No points are present in the sky and on reflective surfaces. Moreover, the point density, which depends on the variable distances between the camera image plane and the positions of the lidar sensor, is generally significantly smaller than the pixel resolution. We can therefore talk about sparse input depth images (see for example Figure 1a, showing the low density of lidar points from the ground).

2. Visibility (hidden parts appear): since points that are not visible from the image view point (hidden points) can be occasionally "seen" by the moving lidar sensor, erroneous values referring to such points can appear in the input depth image. This occurs even when a Z-buffer approach [START_REF] Greene | Hierarchical Z-buffer visibility[END_REF] is used, i.e. only the closest depth values for each pixel are kept (in case multiple values end up in the same pixel location). E.g., Figure 1b shows that depth values from the building behind appear as foreground points.

3. Occlusions (visible parts disappear): for the same reason as above, i.e.

the different acquisition timing and geometry between image and lidar sensors, surfaces normally visible from the image view point do not get a corresponding depth. This can happen when the lidar sensor suffers occlusions at a given instant or because of the scene dynamics. E.g., in Figure 1c, a moving bus that is not present at the moment of the image shot happens to appear in the depth image.

While there is a variety of methods in the literature that deal with the first issue, i.e. that aim at upscaling an irregular input depth image possibly with the guidance of a corresponding color image, little work has been performed

to address the last two issues. In this paper, while inpainting the input depth image, we also intend to tackle the visibility problem. Moreover, we treat at the same time an additional input: a sparse reflectance image derived in the same way as the input depth image (i.e., by naively projecting the lidar point cloud, considering the reflectance information carried out by each point). We will show that the simultaneous use of a reflectance image, which is inpainted jointly with the depth, improves the quality of the produced depth image itself. To jointly inpaint depth and reflectance and concurrently evaluate the visibility of each point (i.e. establish if a single point is reliable or, since non-visible, must be discarded), we formulate an optimization problem with three variables to estimate: depth, reflectance and a visibility attribute per pixel. The inpainting process is also guided by the available color image, by means of a two-fold coupled total variation (TV) regularizer.

The remainder of the paper is organized as follows. In Section 2, we present our approach and mention the related works, in particular on the image-guided depth inpainting problem. In sections 3 and 4 we describe the model used and the primal-dual optimization algorithm that arises, respectively. Finally, in Section 5 we bring experimental evidence that proves the effectiveness of the proposed approach. In the literature there is a variety of methods that aim at upscaling or inpainting an original sparse depth image. Most of them are presented in the context of ToF cameras; thus, a high quality color image is acquired at the same time and can be exploited. We refer to this problem as image-guided depth inpanting. The typical assumption, when exploiting the available image, is that image edges are related to depth edges. Following this principle, many approaches have been proposed, such as methods using different versions of multilateral filtering [START_REF] Chan | A Noise-Aware Filter for Real-Time Depth Upsampling[END_REF][START_REF] Garcia | Pixel weighted average strategy for depth sensor data fusion[END_REF][START_REF] Yang | Fusion of median and bilateral filtering for range image upsampling[END_REF], methods based on Markov Random Fields [START_REF] Diebel | An application of Mmarkov random fields to range sensing[END_REF], and methods using Non-Local Means [START_REF] Huhle | Fusion of range and color images for denoising and resolution enhancement with a non-local filter[END_REF][START_REF] Park | High Quality Depth Map Upsampling for 3D-TOF Cameras[END_REF]. Another family relates to recent methods that make use of optimization [START_REF] Ferstl | Image Guided Depth Usampling using Anisotropic Total Generalized Variation[END_REF][START_REF] Harrison | Image and Sparse Laser Fusion for Dense Scene Reconstruction[END_REF][START_REF] Liu | Guided Depth Enhancement via Anisotropic Diffusion[END_REF][START_REF] Schneider | Semantically Guided Depth Upsampling[END_REF][START_REF] Schwarz | Depth Map Upscaling Through Edge Weighted Optimization[END_REF]. Among these, in [START_REF] Harrison | Image and Sparse Laser Fusion for Dense Scene Reconstruction[END_REF], a method to assign image pixel with a range value, using both image appearance and sparse laser data, is proposed. The problem is posed as an optimization of a cost function encapsulating a spatially varying smoothness cost and measurement compatibility. In the same spirit, the authors of [START_REF] Ferstl | Image Guided Depth Usampling using Anisotropic Total Generalized Variation[END_REF] present an optimization-based depth upsampling method, which uses an Anisotropic Total Generalized Variation (ATGV) term to regularize the solution while exploiting the color image information. Another recent algorithm for the upsampling of sparse depth data is presented in [START_REF] Schneider | Semantically Guided Depth Upsampling[END_REF]. The key idea here is to exploit additional object boundary cues (via structured edge detection and semantic scene labelling) together with usual intenisty cues in a unique optimization framework.

Problem addressed and related work

While presenting good results on images that are not particularly "problematic", in none of the mentioned methods the visibility issue is directly tackled, i.e. there is no explicity estimation of input depth measures to possibly exclude from the inpainting process.We instead intend to estimate visibility, to be able to cope with realistic depth images. To this end, we build on our previous work on lidar-based depth inpainting [START_REF] Bevilacqua | Visibility Estimation and Joint Inpainting of Lidar Depth Maps[END_REF].

W.r.t. the latter, the model is significantly modified to include a reflectance image as well into a new optimization framework. We will show that depth and reflectance mutually benefit of each other in the inpainting process, thus leading to better output results for both. In the next section we present the novel model.

Model

Let Ω ⊆ R 2 be the "full" image support, and Ω S ⊆ Ω the sparse image support where the input images are defined (i.e., there is at least one lidar point ending up there after projection). Given an input depth image v takes values in [0, 1], where v = 0 stands for "hidden" and v = 1 means that the point is visible from the considered image view point. The model considered consists of four terms:

min u∈[um,u M ] r∈[rm,r M ] v∈[0,1] F (u, v|u S ) + G(r, v|r S ) + H(v|u S , r S ) + R(u, r|w) . (1)
F (u, v|u S ) and G(r, v|r S ) are two data-fidelity terms, for depth and reflectance respectively. In both of them the visibility attribute v intervenes.

H(v|u S , r S ) is a term depending exclusively on v, which represents the total cost of classifying input pixels as non-visible. Finally, R(u, r|w) is a regularization term that penalizes the total variation of u and r, by also taking into account the color image w. In the next sections we will detail all the terms composing (1).

Visibility-weighted data-fidelity terms

The data-fitting terms in (1) are meant to enforce fidelity with the original values of depth and reflectance, u S and r S respectively. Deviations from the original values are more penalized if the points are considered "trustful"; conversely, for erroneous original measures (e.g., referring to hidden points) larger deviations are allowed. Therefore we use the visibility attribute v to weight the data terms. For the reflectance data-fidelity term G(r, v|r S ) we have the following expression:

G(r, v|r S ) = η 2 Ω S v|r -r S | dx 1 dx 2 , (2)
where η 2 is a coefficient weighting the term within the model, and dx 1 and dx 2 express the differential lengths in the two image directions. Note that in

(2) an 1 -norm error is used. The 1 norm is considered in substitution of the classical 2 measure of the error for its effectiveness in implicitly removing impulse noise with strong outliers [START_REF] Nikolova | A Variational Approach to Remove Outliers and Impulse Noise[END_REF] and its better contrast preservation [START_REF] Chan | Aspects of Total Variation Regularized L1 Function Approximation[END_REF]. As said, weighting by v relaxes the dependence on the input data for those points classified as hidden.

The depth data-fidelity term, weighted by the coefficient η 1 , is further divided into two terms, as follows:

F (u, v|u S ) = η 1 Ω S max(0, u -u S ) dx 1 dx 2 + Ω S v(max(0, u S -u)) dx 1 dx 2 = F 1 (u|u S ) + F 2 (u, v|u S ) .
(3)

The basic idea behind this separation is to treat differently over-and underestimated depths. Points for which the estimated depth is greater than the original value (u > u S ) most likely correspond to correct input measures, where the over-estimation would be due to the surrounding presence of larger erroneous depths. The expression max(0, u -u S ) is meant to select this kind of points (over-estimated depths). As they are considered reliable, an unweighted data-fitting term, F 1 (u|u S ), is imposed. It is easy to see that for these points the visibility attribute v tends to converge to 1, i.e. they are the best candidates for being classified as visible points. Conversely, the hidden points to remove are sought among depth values which undergo underestimation (u < u S ). These points are taken into account in the second term

F 2 (u, v|u S ),
where the 1 error is weighted by the visibility attribute. Ideally, a fraction of them, the most "problematic" ones, will be classified as hidden (v = 0) and thus not considered in the data fitting cost. Figure 4 shows graphically the depth data-fidelity cost as a function of u -u S . Depending on the value of the visibility attribute v, the 1 -type error |u -u S | is relaxed for negative depth deviations (u < u S ).

Removal cost

The second term of the model ( 1) is meant to penalize the total number of hidden points.

H(v|u S , r S ) = Ω S α(u S , r S )(1 -v) dx 1 dx 2 . ( 4
)
The cost of a single pixel exclusion is proportional to 1 -v, i.e. we have the highest cost for an input pixel when it is totally excluded in the data-fitting cost (v = 0). We individually weight each removal cost, in order to give different importance to each decision visible/hidden. Individual weighting is

given by a coefficient dependent on the original depth and reflectance values,

α(u S , r S ). We generally choose α = k 1 u S + k 2 r S .
The linear dependence of α on the depth and the reflectance "balances" the three terms of (1) depending on v, such that k 1 and k 2 appear to be constants. We will discuss their choice later in this paper (Section 5.1).

Coupled Total Variation

Depth upsampling/inpainting methods that exploit corresponding camera images often relate image edges to depth edges. This has been shown to improve the quality of the reconstructed depth images.

To couple two images in a total variation framework, we adopt the coupled total variation (coupled TV) of [START_REF] Pierre | Luminance-Chrominance Model for Image Colorization[END_REF]:

TV λ (a, b) = Ω (∂ x 1 a) 2 + (∂ x 2 a) 2 + λ 2 (∂ x 1 b) 2 + λ 2 (∂ x 2 b) 2 dx 1 dx 2 . ( 5
)
where λ is a coupling parameter. When λ = 0 the minimization of TV λ encourages the gradient "jumps" to occur at the same locations in a and b.

The coupled TV is then a way to align the edges of an image with those of a given one.

In our problem we have three types of images: a color image w, a depth image u, and a reflectance image r. Figure 3 reports in the bottom row an example of gradient magnitudes related to three images. The gradients of the input depth and reflectance images have been computed after initial interpolation of the latter. As we can clearly see from the image, the color image gradient particularly matches the reflectance one, while being rather dissimilar to the depth gradient. In turn, the reflectance gradient shares some patterns, yet less prominently, with the depth one. See, e.g., the area at the base of the column, where multiple layers mix and produce a similar effect in the two gradient images. We therefore propose to match the three gradients two by two: depth with reflectance, and the same reflectance with the fixed color image. By using the previous definition of coupled TV (5),

we express the regularization term as follows:

R(u, r|w) = TV λ 1 (u, r) + TV λ 2 (r, w) . (6)
After detailing all the terms, our model (1) can therefore be rewritten as follows, the four terms being still distinct:

min u∈[um,u M ] r∈[rm,r M ] v∈[0,1] η 1 Ω S max(0, u -u S ) + Ω S v(max(0, u S -u)) F : Data-fidelity for Depth + η 2 Ω S v|r -r S | G: Data-fidelity for Reflectance + Ω S α(u S , r S ) (1 -v) H: Removal cost + TV λ1 (u, r) + TV λ2 (r, w) R: TV regularization . ( 7
)
In the next section we detail a primal-dual approach to solve (7).

Algorithm

The optimization problem (7) turns out to be convex, but not smooth, due to 1 -type data-fidelity terms, F (u, v|u S ) and G(r, v|r S ), and the total variation regularization term R(u, r|w). Recently, in [START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF] a primal-dual first-order algorithm has been proposed to solve such problems. In Section 4.1 we provide the necessary definitions for the algorithm, which is subsequently described in Section 4.2.

Discrete setting and definitions

Images, considered in Section 3 as continuous functions in R 2 , are here converted into real finite-dimensional vectors. Let M and N be the image dimensions in this discrete setting, and (i, j) the indices denoting all possible discrete locations in the Cartesian grid of size

M × N (1 ≤ i ≤ M , 1 ≤ j ≤ N ).
We then have u, u S , r, r S , v, w, and α ∈ X = R M N , where X is a finite dimensional vector space equipped with a standard scalar product:

u, v X = 1≤i≤M 1≤j≤N u i,j v i,j , u, v ∈ X . (8)
The gradient of an image u ∈ X, ∇u, is a vector in the vector space X 2 with two components per pixel:

(∇u) i,j = ((∇ V u) i,j , (∇ H u) i,j ) . (9)
We compute the gradient components via standard finite differences with Neumann boundary conditions, i.e.:

(∇ V u) i,j = u i+1,j -u i,j i < M 0 i = M (∇ H u) i,j = u i,j+1 -u i,j j < N 0 j = N . ( 10
)
From the definition of gradient, it follows the expression of discrete coupled total variation, which matches the continuous one (5):

TV λ (a, b) = 1≤i≤M 1≤j≤N (∇ H a i,j ) 2 + (∇ V a i,j ) 2 + λ 2 (∇ H b i,j ) 2 + λ 2 (∇ V b i,j ) 2 . (11)
As first suggested by [START_REF] Chan | A Nonlinear Primal-Dual Method for Total Variation-Based Image Restoration[END_REF], a total variation optimization problem can be recast into a primal-dual form that makes its solution easier, by rewriting the gradient norm by means of a vector-valued dual variable. To this end, in our case we first define a "coupled gradient" operator

K λb : X → Y (Y = X 4
), which, applied to an image a ∈ X, expands its gradient to include the one of a reference image b according to a coupling parameter λ. I.e., we have the following element-wise definition:

(K λb a) i,j = ((∇ H a) i,j , (∇ V a) i,j , λ(∇ H b) i,j , λ(∇ V b) i,j ) . ( 12
)
The coupled gradient operator K λb can be further decomposed as

K λb = K + β λ (b)
, according to the following element-wise definition:

(K λb a) i,j = ( Ka) i,j + (β λ (b)) i,j = ((∇ H a) i,j , (∇ V a) i,j , 0, 0) + (0, 0, λ(∇ H b) i,j , λ(∇ V b) i,j ) . ( 13
)
K is the usual gradient operator "padded" with two zero components and it is linear in a; β λ (b) is a bias term, depending on the gradient of the fixed variable b, which determines the last two components of the global coupled gradient operator.

Thanks to the definitions above, we can express alternatively the coupled total variation (11), by introducing the dual variable p ∈ Y :

TV λ (a, b) = max p∈Y K λb a, p Y -δ P (p) = max p∈Y Ka, p Y + β λ (b), p Y -δ P (p) , (14)
where the scalar product in Y is defined as

p, q Y = 1≤i≤M 1≤j≤N
p 1 i,j q 1 i,j + p 2 i,j q 2 i,j + p 3 i,j q 3 i,j + p 4 i,j q 4 i,j p = (p 1 , p 2 , p 3 , p 4 ), q = (q 1 , q 2 , q 3 , q 4 ) ∈ Y , δ P denotes the indicator function of the set

P δ P (p) = 0 if p ∈ P +∞ if p / ∈ P , (15)
and the feasibility set P for the dual variable p, is defined as

P = p ∈ Y | p i,j 2 ≤ 1, ∀i, j , (16)
i.e. p ∞ ≤ 1.

We can now finally express the regularization term of our model R(u, r|w) (6) as the maximization over two dual variables. We then have:

R(u, r|w) = max p∈Y max q∈Y K λ1r u, p Y + K λ2w r, p Y -δ P (p) -δ Q (q) = max p∈Y max q∈Y Ku, p Y + β λ1 (r), p Y + Kr, q Y + β λ2 (w), q Y -δ P (p) -δ Q (q) . ( 17
)
This will let us formulate a discrete version of our joint inpainting problem (7), which falls into the primal-dual optimization framework. As for the other terms in (7), rewritten in discrete notation, we have:

F 1 (u|u S ) = η 1 1≤i≤M 1≤j≤N Φ i,j max(0, u i,j -u S i,j ) F 2 (u, v|u S ) = η 1 1≤i≤M 1≤j≤N Φ i,j v i,j max(0, u S i,j -u i,j ) G(r, v|r S ) = η 2 1≤i≤M 1≤j≤N Φ i,j v i,j |r i,j -r S i,j | H(v|u S , r S ) = 1≤i≤M 1≤j≤N Φ i,j α i,j (1 -v i,j ) (18)
where Φ is a binary mask indicating the initial known pixels, i.e. belonging to the sparse image support Ω S .

A primal-dual algorithm

Thanks to the previous definitions, we can express our model ( 7) in the form of the following saddle-point problem, which is an extension (including two extra variables) of the one presented in [START_REF] Pierre | Luminance-Chrominance Model for Image Colorization[END_REF]:

min u∈X min r∈X min v∈X max p∈Y max q∈Y { K 1 u, p + K 2 r, q -D * 1 (p) -D * 2 (q) +A(u) + B(r) + a(u, v) + b(r, v) + C(v)} . ( 19
)
It is a primal-dual problem with three primal variables (u, r, and v) and two dual variables (p and q) that evolve independently. Each dual variable is particularly linked to the gradient of a primal variable, i.e. p to u, and q to r. D * 1 , D * 2 , A, B, and C are convex functions; a and b are convex w.r.t. each of its respective variables. Globally, the functional is not convex w.r.t. the triplet (u, r, v). By relating ( 7) and ( 19), and using the primal-dual expression of the regularization term reported in (17), we have the following equivalences:

• K 1 u = Ku; • K 2 r = Kr; • D * 1 (p) = -β λ 1 (r), p Y + δ P (p); • D * 2 (q) = -β λ 2 (w), q Y + δ Q (q); • A(u) = F 1 (u|u S ) + δ [um,u M ] (u); • B(r) = δ [rm,r M ] (r); • a(u, v) = F 2 (u, v|u S ); • b(r, v) = G(r, v|r S ); • C(v) = H(v|u S , r S ) + δ [0,1] (v).
An algorithm to solve ( 19) can be derived within the primal-dual optimization framework of [START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF] Algorithm 1 involves the computation of the adjoints to the linear operators K 1 and K 2 (the "zero-padded" gradient operators). It is known that the adjoint of the gradient operator is the negative divergence operator (∇ * = -div). In our case, the adjoint to the operator K 1 : X → Y is a linear operator K * 1 : Y → X consisting in the negative divergence computed only on the two first components of a four-component dual variable p ∈ Y , and by taking finite differences in the opposite direction than the gradient operator (10). These components are in fact the ones related to the primal variable to which the coupled gradient operator has been applied. We then have the following element-wise definition for K * 1 p (the same definition stands for K * 2 q):

(K * 1 p) i,j = -      p 1 i,j -p 1 i-1,j if 1 < i < M p 1 i,j if i = 1 -p 1 i-1,j if i = M -      p 2 i,j -p 2 i,j-1 if 1 < j < N p 2 i,j if j = 1 -p 2 i,j-1 if j = N . (20)
Algorithm 1 Primal-dual based algorithm for depth and reflectance joint inpainting.

1: Inputs: u 0 , r 0 , w, σ, ρ, τ 2: Initialize: u 0 , ū0 ← u 0 , r 0 , r0 ← r 0 , v 0 i,j ← 0.5, p 0 ← (∇u 0 , λ 1 ∇r 0 ), q 0 ← (∇r 0 , λ 2 ∇w) 3: for n = 0, 1, . . . do 4:

p n+1 ← prox σD * 1 (p n + σK 1 ūn ) 5: q n+1 ← prox σD * 2 (q n + σK 2 rn ) 6: v n+1 ← prox ρa(ū n ,•)+ρb(r n ,•)+ρC (v n ) 7: u n+1 ← prox τ A+τ a(•,v n+1 ) (u n -τ K * 1 p n+1 ) 8: r n+1 ← prox τ B+τ b(•,v n+1 ) (r n -τ K * 2 q n+1 ) 9: ūn+1 ← 2u n+1 -u n 10:
rn+1 ← 2r n+1 -r n 11: end for Closed-form expressions for the update rules in Algorithm 1 can be easily computed by applying the definition of proximal operator (see Appendix A). The resulting expressions are reported here below, where P denotes the projection operation over a given real interval, i.e. values are clipped if exceeding the interval limits. Details about the derivation of these expression can be found in Appendix A.

prox σD * 1 (p) = p + σβ λ 1 (r) max(1, p + σβ λ 1 (r) 2 ) (21) prox σD * 2 (q) = q + σβ λ 2 (w) max(1, q + σβ λ 2 (w) 2 ) (22) prox ρa(ū,•)+ρb(r,•)+ρC (ṽ) =      P [0,1] (ṽ) if Φ i,j = 0 P [0,1] (ṽ + ρα -ρη 2 |r -r S |) if Φ i,j = 1, ūi,j ≥ u S i,j P [0,1] (ṽ + ρα -ρη 1 (u S -ū) -ρη 2 |r -r S |) if Φ i,j = 1, ūi,j < u S i,j (23) prox τ A+τ a(•,v) (ũ) =          P [um,u M ] (ũ) if Φ i,j = 0 P [um,u M ] (ũ -τ η 1 ) if Φ i,j = 1, ũi,j > u S i,j + τ η 1 P [um,u M ] (ũ + vτ η 1 ) if Φ i,j = 1, ũi,j < u S i,j -vτ η 1 P [um,u M ] (u S ) otherwise (24) prox τ B+τ b(•,v) (r) =          P [rm,r M ] (r) if Φ i,j = 0 P [rm,r M ] (r -vτ η 2 ) if Φ i,j = 1, ri,j > r S i,j + vτ η 2 P [rm,r M ] (r + vτ η 2 ) if Φ i,j = 1, ri,j < r S i,j -vτ η 2 P [rm,r M ] (r S ) otherwise (25)
The operations indicated in the proximal operators are pixel-wise, although the pixel coordinates have not been made explicit for clearer reading. and comparisons, in Section 5.1 we motivate some critical choices in terms of model and algorithmic parameters.

Experimental results

The algorithm presented in

Parameters of the algorithm and model choices

Our finally resulting joint inpainting model ( 7) consists of four terms:

two data-fidelity terms, F (u, v|u S ) and G(r, v|r S ), a "removal" cost depending solely on the variable v, H(v|u S , r S ), and the two-fold regularization term R(u, r|w). As discussed in Section 3.1, for the data-fidelity terms we opt for an 1 measure of the error, in order to promote more contrasted solutions [START_REF] Chan | Aspects of Total Variation Regularized L1 Function Approximation[END_REF]. The visibility attribute v weights the data matching cost of each single pixel (data matching is more and more relaxed, as v tends to zero, i.e. when that particular point is considered to be excluded). However, over-estimated depths (u > u S ) are not weighted by v but are fully penalized. These values relate to pixels where either there is noise on a visible point that is slightly corrected (u -u S is small), or the value u S represents an outlier (e.g. it is due to a mobile object). At present, we do not have a way to handle the latter case.

In H(v|u S , r S ) (4), each point removal cost is the product between (1 -v)

(the level of "invisibility" of the point) and a coefficient α depending on the local input depth and reflectance: α = k 1 u S + k 2 r S . This choice has been made in order to balance all terms in ( 7) where v appears. Let us now observe the "complete" update rule for v (last case of (23), i.e. for points with under-estimated depth). According to it, we have that at each iteration v is incremented/decremented by a quantity ∆v = ρ (α -η 1 ∆u -η 2 ∆r). Let us suppose that the fluctuations on depth are significantly larger than the fluctuations on reflectance (the appearance of a hidden point can cause a big "jump" in depth, while the reflectance values might still be similar. For the sake of simplicity we can then adjust the value of α only on the basis of the depth input value. The proposed simplified expression for α is then:

α = ku S . (26)
With the assumptions made we therefore have ∆v ∝ (ku S -η 1 ∆u). The attribute v for a certain pixel increases (it gets a higher confidence as a visible point) if ∆u u S < k η 1 , i.e. if the relative depth deviation is below a certain threshold. k is an adimensional parameter that contributes determining this threshold. Conversely, v decreases for relative depth deviations exceeding the threshold. As for the update of v for points with over-estimated depths (second case of ( 23)), if we hypothesize that α, adjusted on depth, is large enough w.r.t. the reflectance deviation, we have that v progressively tends to one (unless large absolute reflectance deviations occur).

As for the regularization term R(u, r|w), we proposed in Section 3.3 to combine two distinct coupled total variation terms: TV λ 1 (u, r) (depth is individually coupled with reflectance) and TV λ 2 (r, w) (reflectance is individually coupled with the color image). By having two separate coupled TV terms, each one encoded by a dual variable that evolves independently from the other one, the reflectance gradient is constantly brought back to the reference gradient of the color image. At the same time the "correct" gradient information is transferred to the depth via the second term. Figure 5 shows an example of results obtained with the algorithm for the same test case as Figure 3. If we observe the input sparse depth image of Figure 3, we see that the major problems come from the fact that depth values referring to the building behind the column appear mixed with foreground depths. With our algorithm we are able to resolve these conflicts, as we can see in the inpainted depth image (Figure 5a). Part of the input pixels have in fact been removed, i.e. classified as non-visible (v = 0). Figure 5c reports the locations of such points in the original depth image. From the histogram of the values of v (Figure 5f) it is evident that the algorithm produces a bipartition of the points according to their visibility attribute. Figure 5 shows also the inpainted reflectance and the final depth and reflectance gradients.

By comparing the latter to the original gradients (Figure 3), we can observe that they end up incorporating elements of the color image gradient, while removing erroneous edges. In the next section we will present more results obtained with our algorithm, also in comparison to other inpainting methods.

Results with urban data

We consider a data set acquired by a MMS system [START_REF] Paparoditis | Stereopolis II: A multi-purpose and multi-sensor 3D mobile mapping system for street visualisation and 3D metrology[END_REF] Our algorithm, presented in Section 4, gives as output the two inpainted images u and r. As for the produced depth image, our algorithm is visually compared with nearest neighbor (NN) interpolation, the anisotropic total generalized variation (ATGV ) method of [START_REF] Ferstl | Image Guided Depth Usampling using Anisotropic Total Generalized Variation[END_REF], and our previous depth inpainting method [START_REF] Bevilacqua | Visibility Estimation and Joint Inpainting of Lidar Depth Maps[END_REF], which does not rely on reflectance information. We refer to the latter as Depth Inpainting with Visibility Estimation (DIVE ). The optimization problem of DIVE is the following:

min u∈[um,u M ] v∈[0,1] η Ω S (max(0, u -y)) 2 dx 1 dx 2 + η Ω S v(max(0, y -u)) 2 dx 1 dx 2 + Ω S (ku S ) 2 (1 -v) dx 1 dx 2 + TV λ (u, w) . ( 27
)
The DIVE problem can be related to our proposed model ( 7), if we consider in the latter η 1 = η, η 2 = 0, λ 1 = λ, and we suppress the coupled TV term related to the reflectance (depth is instead coupled directly with the color image). Moreover, in ( 27) we have a 2 -norm data fidelity term; as a consequence of that, the coefficient of the removal cost term follows a quadratic law (we have α = (ku S ) 2 , instead of α = ku S , as in ( 7)).

As for the produced reflectance image, our algorithm is compared with nearest neighbor (NN) interpolation, the ATGV method of [START_REF] Ferstl | Image Guided Depth Usampling using Anisotropic Total Generalized Variation[END_REF] applied to reflectance, and a reduced version of our model ( 7) limited to reflectance. We refer to this method as Reflectance Inpainting with Visibility Estimation (RIVE ). The RIVE method is derived from the solution of the following optimization problem: r, w) . (28) Also in this case we can derive the considered problem (RIVE) as a simplified version of our proposed model ( 7), where η 1 = 0, η 2 = η, λ 2 = λ, and the coupled TV term related to depth is suppressed. Moreover, the coefficient of the removal cost, while still following a linear law, here depends on the input reflectance r S .

min r∈[rm,r M ] v∈[0,1] η Ω S v|r -r S | dx 1 dx 2 + Ω S (kr S )(1 -v) dx 1 dx 2 + TV λ (
The four examples reported show the better performance of our algorithm in generating complete depth and reflectance images from real lidar measures.

Results with the image Column1, reported in Figure 6, particularly prove the effectiveness of our algorithm in detecting and removing hidden points appearing in the front, thus producing inpainted images correct from the image view point. These points, in yellow/orange according to the color code used for depth, appear mixed to visible points belonging to the column and the fence. By looking at the depth images generated (row (b)), our algorithm is the only one which is able to remove the misleading points and correctly reconstruct the foreground depth plane. This is even more visible by observing the main marble pole highlighted in the 3-D views (row (c)). While other methods are not able to reconstruct the pole, since "distracted" by the interfering background depths, the reconstruction is better performed in our case. Results on the reflectance image confirm the trend. By observing again the main marble pole, we clearly see that the reflectance is better inpainted. This is possible thanks to the joint use of depth information, which helps detecting hidden points by leveraging depth over-and under-estimations, and the coupling with the color image gradient, which helps correctly restoring the edges. Similar considerations can be made for the image Column2 (visual results are reported in Figure 7). Here the box overlaid on the 3-D views indicates an area where points, non-visible from the reference image view point, should be removed. The removal of these points, as well as the inpainting of depth and reflectance, is performed more efficiently by our method.

Figures 8 and9 show results w.r.t. two other images taken peripherally to the scene. For the image Buildings1, we can observe that with our algorithm the inpainted depth and reflectance images look more satisfactory, the pole on the left being completely unveiled as a foreground element. The box overlaid on the 3-D views highlights a part of the scene where the depth values of two trees interfere. Our proposed algorithm (as well as the DIVE method [START_REF] Bevilacqua | Visibility Estimation and Joint Inpainting of Lidar Depth Maps[END_REF]) makes a correct distinction between the two depth layers. Figure 9, reporting results related to the image Buildings2, presents the problem of wrong lidar measures appearing in the front. Our method turns out to be the most effective in clearing out these points, as also shown in the area highlighted by the box.

Performance on visibility estimation

While in the previous section we evaluated the performance of the algorithm in terms of produced inpainted images u and r, we now want to assess the quality of the third output of the algorithm, i.e. v, the visibility attribute.

As visibility is estimated while performing the depth and reflectance estimation, we can say that our algorithm fuses two problems: hidden point removal (HPR) and inpainting. Typically HPR is, instead, possibly performed as a preliminary operation. For HPR "stand-alone" the state of the art is represented by variations of [START_REF] Katz | Direct Visibility of Point Sets[END_REF] that relate the visible point set to the convex hull of a viewpoint-dependent transformation of it, discarding points based on a concavity threshold as seen from the view point.

While this approach is effective, there is in general no globally satisfactory concavity threshold that would both correctly detect hidden surfaces and keep background points close to foreground silhouettes. To compare the two strategies for estimating visibility (the dedicated operation of [START_REF] Katz | Direct Visibility of Point Sets[END_REF] and our "soft" estimation), we show an example in Figure 10, related to the image Column1. In our case, we consider hidden points those depth values that are assigned v = 0 at the end of the algorithm. As for [START_REF] Katz | Direct Visibility of Point Sets[END_REF], a concavity parameter equal to 4 has been chosen after tuning.

The images obtained show that the "quality" of the visibility estimation process is comparable, if not higher with our method. If we observe closely the zoomed-in areas in Figure 10, in fact, we can see that the HPR method wrongly selects points around the silhouettes (see first patch), while sometimes missing the detection of actual hidden points (see last two patches).

As a further test, we also compare our method (which jointly performs visibility estimation and inpainting), with a two-step approach, where visibil- ity estimation (hidden point removal) is performed as a preliminary operation by the algorithm of [START_REF] Katz | Direct Visibility of Point Sets[END_REF]. Depth is subsequently inpainted with the ATGV-based algorithm of [START_REF] Ferstl | Image Guided Depth Usampling using Anisotropic Total Generalized Variation[END_REF]. Figure 11 reports results for such comparison with two images, the two-step approach being denoted as "HPR + ATGV".

In the two cases of Figure 11, we can observe a better outcome with our algorithm. For the image Column1, the preliminary point removal operation is not able to remove all the ambiguities in the central part of the image,

where the depth values of the fence and the column are confused. For the image Buildings2, the HPR method of [START_REF] Katz | Direct Visibility of Point Sets[END_REF] exceeds in removing several points along the upper board of the image, causing blurred edges in the final reconstructed depth image. Besides the benefits observable in the qualitative assessment, the joint approach of our method has the advantage of not requiring an explicit parameter to be globally set (the concavity threshold in the case of [START_REF] Katz | Direct Visibility of Point Sets[END_REF]) to perform HPR. This is instead done in a "soft" way that adapts to the input image.

Quantitative evaluation with a benchmark data set

In this section we perform a quantitative analysis of our algorithm by using the publicly available KITTI Vision Benchmark Suite [START_REF] Geiger | Vision meets robotics: The KITTI dataset[END_REF]. The KITTI suite includes data acquired by a MMS similar to the one considered for our personal data set in Section 5. with two high-resolution color and grayscale video cameras (enabling stereo vision), acquiring images with a pixel resolution equal to 1242 × 375. Accurate 3-D measures are provided by a Velodyne lidar scanner and a GPS localization system. Thus, the lidar measures are generally used as ground truth for algorithm evaluations. In [START_REF] Menze | Object Scene Flow for Autonomous Vehicles[END_REF] a novel data set is presented for stereo benchmarking, which considers also moving objects. By making a special processing on the latter and manually removing erroneous points due to occlusions, ground truth disparity maps are obtained. These maps appear "cleaner" and denser than the input depth images that can be obtained with the raw lidar data, and can therefore be used to evaluate algorithm estimating disparity. To exploit this possibility, as described in (Schneider et al., 2016, Sec. 4.3), we use the ground truth maps of this stereo benchmark data set to have a quantitative evaluation of our depth+reflectance inpainting algorithm. As done by the authors of [START_REF] Schneider | Semantically Guided Depth Upsampling[END_REF], we identify 82 frames (provided ground truth disparity maps) for which we can find correspondences in the raw data set, i.e.

a corresponding color image and related lidar point cloud. We then use the raw data lidar to compute an input depth (e.g., Figure 12a) and use the provided ground truth map to compute a Mean Absolute Error (MAE). The ground truth maps, although denser than the input maps, are sparse, i.e.

they are not defined for all pixels (only about 19% of the pixels have values).

Thus, the MAE is computed only for those pixels which are defined in the respective ground truth map. We computed the MAE for all 82 frames of the found correspondences, for our method and the ATGV-based algorithm of [START_REF] Ferstl | Image Guided Depth Usampling using Anisotropic Total Generalized Variation[END_REF]. As done in Section 5.3, we also compare with a two-step approach, where AGTVbased inpainting is preceded by a hidden point removal (HPR) operation, performed with the algorithm of [START_REF] Katz | Direct Visibility of Point Sets[END_REF]. The resulting average MEAs, which are measured as the average pixel displacement between two disparity maps, are reported in Table 1.

ATGV HPR+ATGV Proposed

Average MAE (px.) 2.13 2.07 1.99 When creating the ground truth maps, the authors of the KITTI benchmark data set have removed objects presenting particular issues in terms of visibility. Other objects are instead manually handled (they are removed from the scene and re-inserted after fitting a CAD model). Thus, the ground truth maps basically consist of the latter and fixed parts of the scene (e.g. streets and walls) that do not yield any ambiguity. Due to this relative "simplicity" of the data set, the performance in terms of average MAE are rather similar among the three methods (ATGV, HPR+ATGV, and proposed method), with our method obtaining a slightly lower error. Nevertheless, we can observe that the ATGV method of [START_REF] Ferstl | Image Guided Depth Usampling using Anisotropic Total Generalized Variation[END_REF] produces more artifacts (see, for example, the reconstructed pole on the left in Figure 12c, in comparison to Figure 12e). Most of these artifacts can be removed by performing a preliminary HPR step (see, in Figure 12b, an example of input depth map cleaned out of ambiguous pixel). The combination of a HPR step and the ATGV-based depth upsampling algorithm of [START_REF] Ferstl | Image Guided Depth Usampling using Anisotropic Total Generalized Variation[END_REF] yields inpainted depth maps with a visual quality comparable to the one of our approach. However, as stated in Section 5.3, with our approach we keep the advantage of having an all-in-one procedure performing jointly inpainting and "soft" visibility estimation (without the need of setting a per-image global threshold as requested by the algorithm of [START_REF] Katz | Direct Visibility of Point Sets[END_REF]). We also expect for our method a greater improvement of the MAE metric and the visual outcome on more complex scenes.

Conclusion

In this paper we presented a novel strategy to jointly inpaint depth and problem is related to disocclusions: the detection of mobile objects is in this case necessary to prevent occlusions in the produced depth and reflectance images. In this case, we could not only use each point measurement as a vote that there is a surface at that point, but we could also exploit the information that the line of sight from the sensor to that point is unobscured, e.g.

as done in [START_REF] Xiao | Street environment change detection from mobile laser scanning point clouds[END_REF]. Other work, concerning modifications to the model presented in this paper, could consist in considering multiple depth candidates per pixel when several points project into the same 2-D location (instead of the Z-buffer approach currently used). = P [0,1] (ṽ) .

(A.4)

! If Φ i,j = 1 and ūi,j ≥ u S i,j , we have:

prox 2 = arg min v 1 2 (v -ṽ) 2 + ρη 2 v|r -r S | + ρα(1 -v) + δ [0,1] (v) = arg min v 1 2 v 2 -vṽ + ρη 2 v|r -r S | -ραv + K + δ [0,1] (v) = arg min v 1 2 v 2 -v (ṽ + ρα -ρη 2 |r -r S |) + K + δ [0,1] (v) = arg min v 1 2 [v -(ṽ + ρα -ρη 2 |r -r S |)] 2 + K + δ [0,1] (v)
= P [0,1] (ṽ + ρα -ρη 2 |r -r S |) .

(A.5) ! If Φ i,j = 1 and ūi,j < u S i,j , we have: if Φ i,j = 0 P [0,1] (ṽ + ρα -ρη 2 |r -r S |) if Φ i,j = 1, ūi,j ≤ u S i,j P [0,1] (ṽ + ρα -ρη 1 (u S -ū) -ρη 2 |r -r S |) if Φ i,j = 1, ūi,j > u S i,j . (A.7) By substituting the optimal value found for u in the splitting condition, we have: u i,j < u S i,j ⇒ ũi,j < u S i,j -vτ η 1 .

! The remaining case is: Φ i,j = 1 and u i,j = u S i,j . This directly implies the solution for the proximal operator: prox 2 = P [um,u M ] (u S ) .

(A.12)

From the previous cases, we can derive the related validity condition on the calculation point ũi,j , i.e.:

-vτ η 1 < ũi,j -u S i,j < τ η 1 .

! Summing up, we have:

prox 3 =          P [um,u M ] (ũ)
if Φ i,j = 0 P [um,u M ] (ũ -τ η 1 ) if Φ i,j = 1, ũi,j > u S i,j + τ η 1 P [um,u M ] (ũ + vτ η 1 ) if Φ i,j = 1, ũi,j < u S i,j -vτ η 1 P [um,u M ] (u S ) otherwise A.15)

Figure 1 :

 1 Figure 1: Examples of parts from a resulting input depth image (bottom row), with the corresponding parts from the reference color image (top row), showing the three issues mentioned: undersampling, appearance of hidden points, and presence of occlusions.

Figure 2

 2 Figure2depicts the scheme of the proposed approach. Given an MMS data set consisting of a lidar point cloud and a set of camera images, we choose among the latter a reference color image (w), and we obtain input depth (u S ) and reflectance (r s ) images by re-projecting the lidar points according to the image geometry. The two lidar-originated images are sparse images with irregular sampling and need to be inpainted. We propose to do that jointly and simultaneously estimate the visibility of the input points,

Figure 2 :

 2 Figure 2: General scheme of the proposed approach. The final outputs of the algorithm are the inpainted reflectance and depth images, r and d respectively, and a binary visibility image v. To represent v, we show the original depth values that finally get v 0.

uFigure 3 :

 3 Figure 3: Example of input depth, reflectance and color images (top row), and their respective gradient images (bottom row). Besides the input depth image, the color map used to encode depth values is reported. Gradients of depth and reflectance are computed on the interpolated versions of the input sparse images, initially obtained by nearest neighbor interpolation.

Figure 4 :

 4 Figure 4: Depth data-fidelity cost F (u, v|u S ) as a function of u -u S for different values of v (η 1 = 1 for simplicity). For over-estimated depths (u -u S > 0) the cost is independent of v, whereas for u -u S < 0 we have different lines as v varies.

  . It consists in a unique loop, where all variables are alternatively updated via proximal operators (see Algorithm 1). The algorithm takes as inputs the initial estimates of the complete depth and reflectance images (u 0 and r 0 , respectively), and the reference intensity image w. It also requires three parameters inherent to the algorithm: σ and τ , which are related to each other by the relation 16τ σ ≤ 1[START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF], and ρ, which is a parameter regulating the update speed of v.

  Section 4 is evaluated with a new data set acquired in an urban scenario by a Mobile Mapping System (MMS), composed of lidar measures and camera-originated images. With this data set, we provide a qualitative evaluation of our algorithm in comparison with other methods, by showing the reconstructed depth and reflectance images, and we assess the quality of the visibility estimation task, which is a crucial characteristic of our algorithm. Moreover, we also provide a quantitative analysis by means of a benchmark data set publicly available. Before showing results

Figure 5 :

 5 Figure 5: Output of the proposed algorithm for the image Column1 : (a) Inpainted depth, (b) Inpainted reflectance, (c) Removed points (v = 0), (d) Final depth gradient, (e) Final reflectance gradient, (f) Final histogram of v.

  Figures 6-9 present results for four images (cropped w.r.t. the full size) of the data set: Column1, Column2, Buildings1, Buildings2. For each reference image, the input sparse depth and reflectance images, obtained via projection, are shown, as well as the inpainted depth and reflectance images, obtained with four different methods. For the output depth images of Figure 8 and 9 we added some shading by modulating the color intensity of each pixel based on the zenith angle of the normal vector, to emphasize high-frequency changes. Moreover, for the inpainted depths, an alternative view of the resulting 3-D point cloud is proposed, where the coordinates of the points are retrieved thanks to the computed depths and color texture is applied to enrich the points. A color box is overlaid to the first of these 3-D views to highlight areas where the comparison between the different methods is particularly significant.

Figure 6 :

 6 Figure 6: Visual results for the image Column1. Row (a) shows the related input images: depth (with a 3-D zoom), reflectance, and reference color image. Rows (b) and (c) report the results obtained in terms of inpainted depth images (with related 3-D zoomed-in view) with the algorithms indicated below. Row (d) shows the inpainted reflectance images obtained with different methods, our proposed method always reported as last.

Figure 7 :

 7 Figure 7: Visual results for the image Column2. Row (a) shows the related input images: depth (with a 3-D zoom), reflectance, and reference color image. Rows (b) and (c) report the results obtained in terms of inpainted depth images (with related 3-D zoomed-in view) with the algorithms indicated below. Row (d) shows the inpainted reflectance images obtained with different methods, our proposed method always reported as last.

Figure 8 :

 8 Figure 8: Visual results for the image Buildings1. Row (a) shows the related input images: depth (with a 3-D zoom), reflectance, and reference color image. Rows (b) and (c) report the results obtained in terms of inpainted depth images (with related 3-D zoomed-in view) with the algorithms indicated below. Row (d) shows the inpainted reflectance images obtained with different methods, our proposed method always reported as last.

Figure 9 :

 9 Figure 9: Visual results for the image Buildings2. Row (a) shows the related input images: depth (with a 3-D zoom), reflectance, and reference color image. Rows (b) and (c) report the results obtained in terms of inpainted depth images (with related 3-D zoomed-in view) with the algorithms indicated below. Row (d) shows the inpainted reflectance images obtained with different methods, our proposed method always reported as last.

Figure 10 :

 10 Figure 10: Detected hidden points in the case of the image Column1, by the state-of-theart method of[START_REF] Katz | Direct Visibility of Point Sets[END_REF] and our method. The three patches below each image represent zoomed-in areas of the images themselves at same locations.

Figure 11 :

 11 Figure 11: Comparison between our joint approach and a two-step approach, where visibility estimation and inpainting are performed separately, on the images Column1 (a) and Buldings2 (b).

Figure 12 :

 12 Figure 12: Case example from the 2015 KITTI stereo benchmark data set. For each input depth map (a), we have a ground truth disparity map available yet sparse (f), w.r.t. which it is possible to compute an error by only considering the pixels where it is defined. By applying a hidden point removal (HPR) algorithm to the input depth data it is possible to create a new input map where background hidden pixels have been removed (b). Results for different depth inpainting strategies are reported (c, d, e).

  reflectance images with the guidance of a co-registered color image, and by simultaneously estimating a visibility attribute for each pixel. The problem studied and the proposed approach are particularly suited for data sets acquired by Mobile Mapping Systems (MMS): vehicles that can easily image urban scenes by means of optical cameras and lidar sensors. By projecting the 3D lidar points onto a chosen reference image, we obtain depth and reflectance images, which suffer of practical issues due to the big diversity of the lidar and optical sensor acquisitions. By estimating visibility, we aim at solving one of these issues, i.e. the appearance (in depth and reflectance) of parts of objects non-visible from the image view point, but captured by the lidar sensor. Those points are meant to be detected by our algorithm and thus discarded in the inpainting process. The proposed approach consists in a variational optimization problem, where three variables (depth, reflectance, and visibility) are simultaneously estimated. As a regularization term, a twofold coupled total variation (TV) term is proposed, where the gradients of depth, reflectance and color image are matched two by two, by leveraging the inherent correlation between them. The proposed algorithm is compared, in terms of inpainted images, to other inpainting algorithms, which do not take into account the simultaneous detection of possibly erroneous measures.The clear superiority of the proposed method w.r.t. the latter proves that the visibility estimation is a necessary step. Another comparison is made with a simplified version of the algorithm, which accounts for visibility but considers alternatively either depth or reflectance. The worse performance of the simplified algorithm indicates that the joint exploitation of depth and reflectance is a key aspect for the success of the algorithm. The mutual benefit comes from the fact that depth is particularly important for the visibility estimation task; in turn, reflectance is crucial in restoring the correct edges, via coupling with the color image. Future work will continue in the direction of solving practical issues with lidar-based images to inpaint. Notably, another

  ṽ) 2 + ρη 1 Φv max(0, u S -ū) + ρη 2 Φv|r -r S | + ραΦ(1 -v) + δ [0,1] (v) (A.3)! If Φ i,j = 0 (point out of the sparse domain), we trivially have: ṽ) 2 + δ [0,1] (v)

  ṽ) 2 + ρη 1 v(u S -ū) + ρη 2 v|r -r S | +ρα(1 -v) + δ [0,1] vṽ + ρη 1 v(u S -ū) + ρη 2 v|r -r S | -ραv +K + δ [0,1] v (ṽ + ρα -ρη 1 (u S -ū) -ρη 2 |r -r S |) (ṽ + ρα -ρη 1 (u S -ū) -ρη 2 |r -r S |)] 2 +K + δ [0,1] (v) = P [0,1] (ṽ + ρα -ρη 1 (u S -ū) -ρη 2 |r -r S |) .

•

  prox 4 = prox τ B+τ b(•,v) r) 2 + τ η 2 Φv|r -r S | + δ [r m ,r M ] (r) (A.14) ! If Φ i,j = 0, we trivially have: r) 2 + δ [rm,r M ] (r) = P [rm,r M ] (r) .

(

Table 1 :

 1 Average Mean Absolute Error (MAE), i.e. average pixel displacement between ground truth and reconstructed disparity maps, obtained by averaging the results of 82 frames of the 2015 KITTI stereo benchmark data set.

$ This study has been carried out with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of the Investments for the future Programme IdEx Bordeaux (ANR-10-IDEX-03-02). J.-F. Aujol also acknowledges the support of the Institut Universitaire de France.

Appendix A. Derivation of the proximal operators in Algorithm 1

In this section we detail the derivation of the closed-form expressions of the proximal operators for the update of three primal variables (v, u, and r) in Algorithm 1, as listed in Section 4.2. Let f : R n → R ∪ {+∞} be a closed proper convex function. The proximal operator or mapping prox f : R n → R n of f [START_REF] Parikh | Proximal Algorithms[END_REF] is defined by:

Broadly speaking, the proximal operator of a function is a mathematical tool that allows to make an approximation to a certain value, while making a compromise between the accuracy of the approximation and a cost given by the function itself.

Given the general definition (A.1), we can derive the expressions for the proximal operators of the functions considered in our algorithm. We have that the operations involved are independent for each coordinate of the processed images. Therefore, the expressions reported below are to be intended per coordinate, although the spatial indices indicating a particular pixel location are not specified for brevity.

(The variable p, even considered for a single coordinate, is vectorvalued.)

• prox 2 = prox ρa(ū,•)+ρb(r,•)+ρC (ṽ)

! If Φ i,j = 0, we trivially have:

(A.9)

! If Φ i,j = 1 and u i,j > u S i,j , we have:

(A.10) By substituting the optimal value found for u in the splitting condition, we have:

! If Φ i,j = 1 and u i,j < u S i,j , we have:

! If Φ i,j = 1 and r i,j > r S i,j , we have:

By substituting the optimal value found for r in the splitting condition, we have:

! If Φ i,j = 1 and r i,j < r S i,j , we have:

(A.17)

By substituting the optimal value found for r in the splitting condition, we have:

! The remaining case is: Φ i,j = 1 and r i,j = r S i,j . This directly implies the solution for the proximal operator:

From the previous cases, we can derive the related validity condition on the calculation point ri,j , i.e.:

|r i,j -r S i,j | < vτ η 2 .

! Summing up, we have:

if Φ i,j = 0 P [rm,r M ] (r -vτ η 2 ) if Φ i,j = 1, ri,j > r S i,j + vτ η 2 P [rm,r M ] (r + vτ η 2 ) if Φ i,j = 1, ri,j < r S i,j -vτ η 2 P [rm,r M ] (r S ) otherwise (A.19)