Elian Dib

Mikaël Le Pendu

Xiaoran Jiang

Christine Guillemot
email: christine.guillemot@inria.frlependum@scss.tcd.ie

Mikael Le Pendu

Super-rays based low

come

Introduction

Light field imaging has emerged as a promising technology for a number of applications going from photo-realistic image-based rendering to a variety of computer vision applications such as e.g. 3D modeling, object detection, classification, recognition. However, light fields represent very large volumes of high dimensional data, hence the need for designing efficient compression algorithms. In this paper, we focus on the problem of compression of dense light fields, as those captured by plenoptic cameras, which represent very large volumes of highly redundant data. While a number of methods have already been published in the literature aiming at adapting standardized solutions (in particular HEVC) to light field data as in [START_REF] Conti | HEVC-based light field image coding with bipredicted self-similarity compensation[END_REF][START_REF] Monteiro | Light field HEVC-based image coding using locally linear embedding and self-similarity compensated prediction[END_REF][START_REF] Li | Efficient intra prediction scheme for light field image compression[END_REF], here we focus on the problem of reducing the angular dimension of light fields with a low rank approximation method. An homography-based low-rank approximation method called HLRA has been shown to give very good light fields compression performances in [START_REF] Jiang | Light field compression with homography-based low-rank approximation[END_REF]. However, the low rank approximation in [START_REF] Jiang | Light field compression with homography-based low-rank approximation[END_REF] is done globally on the entire light field.

In this paper, we explore the use of super-ray based local low rank approximation models for further improving the performance of the light field compression algorithm. The super-rays are used here as a way to better expose redundancy across the different views compared to a global homography-based alignment as done in [START_REF] Jiang | Light field compression with homography-based low-rank approximation[END_REF], and in turn to further reduce the angular light field dimension. The concept of super-ray has been initially introduced in [START_REF] Hog | Super-rays for efficient light field processing[END_REF] as an extension of super-pixels to address the computational complexity issue in light field image processing tasks. Super-rays can be seen as the clustering of rays of the light field within and across views, hence corresponding to the same set of 3D points of the imaged scene. In this paper, a dedicated super-ray construction method is first proposed so that super-rays can be used as the support of the low rank approximation. The central view is taken as a reference view on which super-pixels are constructed using the SLIC algorithm [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF]. The super-rays which are constructed must group super-pixels which are consistent across the views while being constrained to be of same shape and size. The corresponding super-pixels in all views are found thanks to disparity compensation. In order to do so, a novel method is proposed to estimate the disparity for each super-ray using a low rank prior, so that the super-rays are constructed to yield the lowest approximation error for a given rank.

Due to occlusions, individual super-pixels may contain less pixels than the complete area captured in the corresponding super-ray (gathering information from superpixels of all the views). The super-pixels are thus extended with the occluding pixels to satisfy the shape and size constraint. A low rank approximation is then computed for the set of extended super-rays that are stacked in a vectorized form in different columns of a matrix X. The rank r approximation of the set of extended super-rays is expressed as a product of a matrix B, containing r columns corresponding to basis (or eigen) images in a vectorized form, with a matrix C containing weighting coefficients. The matrix B is re-arranged into two sets of eigen images, respectively corresponding to the light rays visible and occluded in the central view. An image sequence formed with both those sets of 'occluded' and 'visible' eigen images is encoded using HEVC.

The method is also compared against one HEVC-based scheme which first forms a pseudo video sequence from all the views and then code this sequence with HEVCinter coding [START_REF] Rizkallah | Impact of light field compression on focus stack and extended focus images[END_REF], and with the verification model of JPEG-Pleno (VM 1.1) [START_REF]Jpeg pleno light field coding vm 1.1[END_REF]. Experimental results show that a rate saving of up to -58.19% and -54.33% can be achieved compared with [START_REF] Rizkallah | Impact of light field compression on focus stack and extended focus images[END_REF] and [START_REF]Jpeg pleno light field coding vm 1.1[END_REF]. This section does not give an exhaustive review of all light field compression techniques proposed so far but rather focuses on the most recent work on the topic. Existing solutions can be broadly classified into two categories: approaches directly compressing the lenslet images or approaches coding the views extracted from the raw data. Methods proposed for compressing the lenslet images mostly extend HEVC intra coding modes by adding new prediction modes to exploit similarity between lenslet images [START_REF] Conti | HEVC-based light field image coding with bipredicted self-similarity compensation[END_REF], [START_REF] Monteiro | Light field HEVC-based image coding using locally linear embedding and self-similarity compensated prediction[END_REF] . In [START_REF] Tabus | Lossy compression of lenslet images from plenoptic cameras combining sparse predictive coding and jpeg 2000[END_REF], a different approach is proposed using a depth based segmentation of the light field into 4D spatio angular blocks with prediction followed by JPEG2000.

Solutions extracting and coding views, code the views as pseudo video sequences using HEVC [START_REF] Rizkallah | Impact of light field compression on focus stack and extended focus images[END_REF], [START_REF] Liu | Pseudo-sequence-based light field image compression[END_REF], or the latest JEM coder [START_REF] Jia | Optimized inter-view prediction based light field image compression with adaptive reconstruction[END_REF], or extend HEVC to multi-view coding [START_REF] Ahmad | Interpreting plenoptic images as multiview sequences for improved compression[END_REF]. Methods using view synthesis have also been proposed in [START_REF] Jiang | Light fields compression using depth image based view synthesis[END_REF] using convolutional neural networks to synthesize all the views from the corner views, or in [START_REF] Zhao | Light field image coding via linear approximation prior[END_REF] where the authors use a linear approximation computed with Matching Pursuit for disparity based view prediction.

The most closely related work to the proposed method is the approach proposed in [START_REF] Verhack | Steered mixture-of-experts for light field coding, depth estimation, and processing[END_REF], in the sense that local models exploiting disparity are considered in the 4D rays space for compression. However, the models in [START_REF] Verhack | Steered mixture-of-experts for light field coding, depth estimation, and processing[END_REF] are 4D steered Gaussian Mixture Models (GMM), referred to as steered mixture of experts, while the proposed local models are based on low rank priors. The second most related work is the method in [START_REF] Jiang | Light field compression with homography-based low-rank approximation[END_REF] introducing a global homography-based low rank approximation approach for light field compression. However, the proposed approach differs from the above approach, and extends it, by proposing a coding scheme based on local low rank models using a novel disparity estimator also exploiting low rank priors.

Notations and scheme overview

Let L(u, v, x, y) be a light ray of the light field L, and (u, v, x, y) its coordinates using the two plane parametrisation, where (u, v) and (x, y) are the angular (view) and spatial (pixel) coordinates respectively. A view with angular coordinates (u, v) is denoted I u,v . The overall scheme is depicted in Fig. 1 and comprises the following steps:

• Super-rays construction with the help of the disparity maps of all views.

• Low-rank based disparity optimization and compensation of each super-ray.

• Super-rays extension so that all super-pixels forming a SR are of same size.

• Low rank approximation of the set of extended super-rays.

• Splitting the eigen images into a first set {B l v } l∈ 1,r (with r the approximation rank) corresponding to the light rays visible in the central view, and a second set {B l o }, very sparse, and containing eigen values of super-rays extensions. • Encoding the information to be transmitted to the decoder, i.e., the two sets of eigen images {B l v } and {B l o }, the matrix of coefficients C, and side information needed at the decoder to find the super-rays and perform the inverse alignment. These different steps are detailed below.

Here, we present a super-ray construction method which constrains the projections of a super-ray on each view to be of same shape and size.

Light field over-segmentation

A k-means clustering method is applied on the entire light field in order to group pixels of all the views which are in local regions and which are similar in colour (in CIELAB colourspace) and in disparity. This clustering, extending the concept of super-pixels, forms so-called super-rays since it groups pixels (and super-pixels) across views which are assumed to correspond to the same set of 3D points in space. This leads to a segmentation of all the views into super-pixels which are coherent across the views. The clustering into super-rays is performed using the disparity estimation method of [START_REF] Jiang | Depth estimation with occlusion handling from a sparse set of light field views[END_REF]. We then assume that all pixels belonging to a given super-ray have the same disparities referred to as the centroid disparities.

To compute the low rank approximation, the super-pixels in each super-ray need to be aligned by disparity compensation. We show in the experimental section that disparity estimated independently of the low rank constraint, e.g. as in [START_REF] Jiang | Depth estimation with occlusion handling from a sparse set of light field views[END_REF], is not optimal for our low rank based coding scheme. We hence propose in the next section a novel disparity optimization method for each super-ray based on a low rank prior.

Disparity estimation of super-rays using a low rank prior

We use a low rank prior on each super-ray as a means to locally approximate the light field and jointly perform disparity compensation so as to align similar pixels from different views and render the approximation even more efficient. For a given super-ray, let h d be the set of disparity compensation operators h d i for each super-pixel at a view (u i , v i ) forming the super-ray. The disparity compensation h d i is defined as

h d i : (x, y) → (x+d x (u i -u c ), y+d y (v i -v c ))
, and d is the vector d x d y to estimate. The disparity estimation algorithm proceeds as follows. The horizontal and vertical disparity values d x and d y are first initialized to "0". Note that even after disparity compensation, the aligned super-pixels forming a super-ray may have different shapes and sizes, in particular due to occlusions, which would make the low rank approximation impractical. To cope with this difficulty, super-pixels with missing information are filled using occluding pixels. The set of aligned and completed super-pixels is called extended super-ray. All super-pixels forming a given extended super-ray SR are stored, in a vectorized form, in a matrix

SR = [SP 1 | SP 2 | ...] , (1)
where SP i is the vectorized version of the super-pixel from view (u i , v i ).

For each extended super-ray, we then solve

argmin d,M SR • h d -M 2 F s.t. rank(M) = r, (2)
where M is a rank r matrix and . F is the Frobenius norm. The notation "•" stands for the application of the disparity compensations to align all the super-pixels forming a super-ray.

We iteratively solve the problem by successively finding M while fixing d, and updating d while fixing M, at each iteration, until convergence. For a fixed d, the matrix M of rank r is obtained as M = UΣ r V T where UΣV T is the singular value decomposition (SVD) of SR • h d and Σ r contains only the r largest singular values of Σ. The set disparity d is updated, by solving (2) for the rank r matrix M fixed. Minimizing (2) is not trivial due to the non linearity of the term SR • h d .

Assuming the change ∆d at each iteration to be small, we can approximate this update by local linearity as

SP i • h d+∆d i ≈ SP i • h d i + ∂[SP i • h d i ] ∂d x ∆d x + ∂[SP i • h d i ] ∂d y ∆d y (3)
= SP i • h d i + ∂[SP i • h d i ] ∂x (u i -u c )∆d x + ∂[SP i • h d i ] ∂y (v i -v c )∆d y (4)
Given this approximation, d is updated with the increment ∆d obtained by solving

∆d = argmin ∆d i R i -J i ∆d 2 F , (5)
with,

R i = M i -SP i • h d i and J i = (u i -u c ) ∂[SP i •h d i ] ∂x (v i -v c ) ∂[SP i •h d i ] ∂y
. This problem has the following analytical solution:

∆d = i J i J i -1 i J i R i (6)
Equation 2 minimizes the approximation error of the matrix formed by the aligned views. This approximation error corresponds to the PSNR that we call PSNRin in the experimental section. We however show in the sequel that significantly better results are obtained by taking into account the inverse warping (i.e. the corresponding interpolation errors) in the minimization. At each iteration, for a given d we evaluate SR -

(h d ) -1 • M 2 F
to eventually keep the disparity values d minimizing the superray reconstruction error, corresponding to a PSNR that we call PSNRout.

Low Rank Approximation

All the matrices corresponding to the different extended and aligned super-rays SR k are stacked in a matrix X of dimension R m×n , where n is the number of views and m is the number of pixels per view of the concatenated extended super-rays. In other words, the matrix X contains all the pixels of all the views that are locally aligned super-ray per super-ray.

The matrix X is factorized into a product of a low-rank matrix B ∈ R m×r and a coefficient matrix C ∈ R r×n , with r ≤ n, as argmin B,C X -BC 2 F , (7) where . F is the Frobenius norm. Optimal factorization is obtained from the singular value decomposition (SVD) of X into U ΣV . Assuming the singular values in Σ are in decreasing order, we take B as the r first columns of U Σ, and C as the r first rows of V . The set of super-rays {SR} forming the entire light field can thus be approximated by a linear combination of columns of B, therefore significantly reducing the amount of data at the cost however of an approximation error.

Compression Scheme

The above low rank approximation leads to a matrix B containing in its columns eigen images of the extended views (set of extended super-pixels) and to a matrix C containing coefficients. For the purpose of encoding the matrix B using HEVC, the data in each column l is re-arranged into two eigen images B l v and B l o . B l v contains the entries derived from the portions of the extended super-rays which are visible in the central view. B l o contains the entries corresponding to areas occluded in the central view. Two segmentation maps M v and M o are also computed so that each region in M v (respectively M o ) is labelled with the index of the super-ray that will be reconstructed from the collocated region in the B l v images (respectively B l o ). Examples of eigen images B l v and B l o with their associated maps M v and M o are illustrated in Fig. 2.

The data in the images B l v and B l o and the matrix C are quantized on 16 bits. The quantized entries of the C matrix are directly transmitted (e.g. coded using a fixed length code) since the corresponding cost is quite negligible. As shown in Fig. 2, the first 4 eigen images B l v and B l o , show some correlation. Hence to code them efficiently using HEVC-Inter, the empty regions of the B l o images are completed with the collocated regions of the corresponding B l v images. The B l v and B l o images are then interleaved as

B 1 v , B 1 o , ..., B 4 v , B 4
o , and coded using HEVC-Inter with a GOP size of 2. The remaining columns of these matrices being less correlated are coded using HEVC-Intra.

It is also necessary to transmit as side information the two centroid disparity values per super-ray as well as the two segmentation maps M v and M o . The two segmentation maps are compressed without loss using the compresso method described in [START_REF] Matejek | Compresso: Efficient compression of segmentation data for connectomics[END_REF]. The disparity information (two values d x , d y per super-ray) is encoded using a fixed length code on 32 bits, which is quite negligible given the small number of values to be transmitted (in the experiments we considered a partition of the light fields into 200 super-rays for natural light fields and 240 for the synthetic ones).

Light field reconstruction

On the receiver side, the decoded sequence of images Bl v and Bl o need to be merged to form the decoded low rank matrix B which contain the eigen images of all the aligned and extended super-rays. This matrix B is multiplied by the matrix C to obtain the matrix X of aligned and extended super-rays. An inverse disparity compensation is then performed per super-ray using its centroid disparities. Each view of the light field is then reconstructed by progressively mapping each super-ray on each view starting from the super-rays having the smallest disparity, i.e. starting from background pixels. Pixels having a higher disparity over-write previously rendered pixels at the same position with a lower disparity.

Depending on the complexity of the scene, on the number of objects and depth layers, some parts in the light field may be occluded by several objects. The corresponding pixels are referred to here as multiply-occluded pixels. One sets of eigen images {B l o } l is not sufficient to represent multiply-occluded pixels. One alternative would be to transmit additional sets {B l o } l to represent these pixels. However, the number of such pixels is quite limited and the corresponding matrix would be very sparse. Instead, we inpaint the corresponding pixels in the reconstructed light field L using a low rank completion matrix. The low rank matrix completion problem is posed as the search for the minimum nuclear norm matrix L with entries equal to those of the matrix L for the known elements of L. The problem is mathematically formulated as min

L || L|| * s.t. ∀(i, j) ∈ Ω, Lij = Lij , (8)
where Ω is the set of indices of the known elements in L, and ||.|| * is the nuclear norm (convex approximation of the rank). This minimization is solved using the Inexact ALM (IALM) technique [START_REF] Lin | The augmented Lagrange multiplier method for exact recovery of corrupted low-rank matrices[END_REF].

Experimental Results

The performances of the proposed disparity estimation and light field compression methods have been evaluated for the luminance component of light fields shown in Fig. 4 and coming from : 1/-INRIA dataset [START_REF]INRIA Lytro image dataset[END_REF] which contains LFs captured by a second generation Lytro Illum camera from which we use the 9 × 9 central subaperture images of 625 × 434 pixels each; 2/-the HCI dataset [START_REF] Wanner | Datasets and benchmarks for densely sampled 4D light fields[END_REF] synthetic LFs 9 × 9 views of 768×768 pixels; 3/-the ICME 2016 Grand Challenge dataset [START_REF]ICME 2016 Grand Challenge dataset[END_REF] from which we take the 9 × 9 central sub-aperture images. The Lytro LFs have been decoded using the Matlab Light Field Toolbox v0.4 [START_REF] Dansereau | Light Field Toolbox for Matlab[END_REF] with gamma correction. The methods involving a low rank prior (HLRA and our method Local Low Rank Approximation (LLRA)) were evaluated for varying target rank values r = {1, 3, 5, 15, 30, 60} and HEVC quality parameter QP ∈ {2, 6, 14, 20, 26, 38}. We retain the (r,QP ) pairs corresponding to the points on the convex enveloppe of the rate distortion graph. For HEVC Lozenge QP = {10, 14, 17, 20, 23, 26}. The disparity maps used for LLRA and JPEG Pleno were generated using [START_REF] Jiang | Depth estimation with occlusion handling from a sparse set of light field views[END_REF]. For JPEG Pleno we adapted the existing F ountain V incent 2 configuration files for use on 9 × 9 LFs and used it for every LF.

Performance analysis of the disparity estimation using a low rank prior

We first assess in Fig. 3 the performance of the proposed super-ray disparity estimation method based on a low rank prior against an exhaustive search, and disparity computed independently of the low rank constraint (in the experiments we used disparity values computed with [START_REF] Jiang | Depth estimation with occlusion handling from a sparse set of light field views[END_REF]). In the exhaustive search approach, the optimal disparity values for each super-ray are obtained by evaluating the low rank approximation error for a wide range of integer, and then fraction-pixel refined, disparity values. We observe that using a low rank prior improves significantly the performance compared to using input disparity maps computed with no rank constraint. Moreover, our method compares favorably with the exhaustive search. We also observe that aiming at maximizing the PSNR after inverse warping (PSNRout) during the SR disparity optimization rather than maximizing the PSNR of the low rank approximation of the aligned SR (PSNRin) allows us to mitigate the information loss due to bi-cubic interpolation in forward and inverse warping. In the sequel, we use the low-rank based and exhaustive search disparity estimation with maximization of outter PSNR.

Performance analysis of the complete scheme

Tab. 1 shows the Bjontegaard rate savings obtained for the HLRA method [START_REF] Jiang | Light field compression with homography-based low-rank approximation[END_REF] and two variants of our method with respect to the HEVC-Lozenge [START_REF] Rizkallah | Impact of light field compression on focus stack and extended focus images[END_REF] and JPEG Pleno Verification Model [START_REF]Jpeg pleno light field coding vm 1.1[END_REF] anchors and we consistently observe rate savings on natural light fields compared with the two anchors. In particular, for real light fields, we observe comparable performance by using the low rank prior to compute the centroid disparities or by performing an exhaustive search. However, for synthetic light fields, the lower performance is due to the initialization of the gradient descent ((d x , d y ) = (0, 0)) that is not appropriate for larger baselines. By initializing the centroid disparities of each super-ray to the median disparity of the pixels of a disparity map [START_REF] Jiang | Depth estimation with occlusion handling from a sparse set of light field views[END_REF], we observe in Fig. 5 a better performance for LLRA for light fields with larger baselines such as StillLife and Buddha. For these light fields, LLRA can outperform HLRA despite the extra signalling cost due to the segmentation in super-rays.

Conclusion

This paper has presented a compression scheme for light fields using super-ray based local low rank models. A novel method for disparity estimation and compensation was proposed so that the super-rays are constructed to yield the lowest approximation error for a given rank. The proposed representation is actually based on two low rank models, one for the central view pixels that are visible in all views and one for the occlusions. We could observe consistent bitrate savings on all natural light fields compared to reference methods (HEVC and JPEG-Pleno). However, compared to global homographies (HLRA), the segmentation in super-rays implies sending extra side information that can be penalizing at low rate, for light fields with small baselines. However, for larger baselines, the performances benefit from the better alignment and approximation of the local model.

Figure 1 :

 1 Figure 1: Overview of the super-ray based low rank approximation and coding chain (encoder), and Light field reconstruction steps (decoder).

Figure 2 :

 2 Figure 2: First 4 eigen-images (shown in false colors) and segmentation maps for the visible (B l v , M v ) and occluded (B l o , M o ) sets of pixels ("butterfly" synthetic light field [17]).

Figure 3 :

 3 Figure 3: PSNR-Rate performance of the low rank based disparity estimation compared with an exhaustive search and the use of a disparity map computed with no rank constraint.

Figure 4 :

 4 Figure 4: Test light fields. Top row (first two): HCI (Buddha, StillLife) [17]. Top row (last three): INRIA Dataset (Bench, Fruits, Toys) [20]; Bottom row: ICME Dataset (Fountain Vincent 2, Friends 1, StonePillars, Vespa) [21].

Figure 5 :

 5 Figure 5: PSNR-Rate performance for two synthetic light fields.

Table 1 :

 1 93% -51.58% -14.41% -8.28% -11.91% -5.65% Friends 1 -65.62% -73.52% -57.44% -43.07% -58.19% -42.87% S P I -67.84% -72.45% -54.95% -54.33% -50.24% -50.73% Vespa -49.57% -61.21% -37.85% -12.61% -36.54% -5.09% Bench -67.37% -64.25% -46.56% -25.56% -45.65% -25.16% Bjontegaard rate savings for HLRA and LLRA compression schemes with respect to HEVC-Lozenge and JPEG Pleno VM anchors.

	Method	HLRA	LLRA (exhaustive	LLRA (low rank
				disparity search)	prior disparity)
	Anchor HEVC	JPEG	HEVC	JPEG	HEVC	JPEG
		Lozenge	Pleno	Lozenge	Pleno	Lozenge	Pleno
	F V 2	-40.

Fruits -55.15% -59.19% -26.43% -12.54% -27.05% -12.74% Toys -58.98% -81.05% -41.32% -48.96% -42.09% -48.90%

This work has been supported by the EU H2020 Research and Innovation Programme under grant agreement No 694122 (ERC advanced grant CLIM).


Mira Rizkallah
email: mira.rizkallah@inria.fr

Thomas Maugey

Christine Guillemot

Graph-based Spatio-angular Prediction for Quasi-Lossless Compression of Light Fields

Graph-based transforms have been shown to be powerful tools for image compression. However, the computation of the basis functions becomes rapidly untractable when the support increases, i.e. when the dimension of the data is high as in the case of light fields. Local transforms with limited supports have been investigated to cope with this difficulty. Nevertheless, the locality of the support may not allow us to fully exploit long term dependencies in the signal. In this paper, we describe a graph-based prediction solution that allows taking advantage of intra prediction mechanisms as well as of the good energy compaction properties of the graph transform. The approach relies on a separable spatio-angular transform and derives low frequency spatio-angular coefficients from one single compressed reference view and from the high angular frequency coefficients. In the tests, we used HEVC-Intra, with QP=0, to encode the reference frame with high quality. The high angular frequency coefficients containing very little energy are coded using a simple entropy coder. The approach is shown to be very efficient in a context of high quality quasi-lossless compression of light fields.

Introduction

The availability of commercial light field cameras has given momentum to the development of light field compression algorithms. Many solutions proposed so far adapt standardized image and video compression solutions (in particular HEVC) to light field data (e.g. [START_REF] Liu | Pseudosequence-based light field image compression[END_REF] [2] [START_REF] Li | Compression of unfocused plenoptic images using a displacement intra prediction[END_REF]). Recent work investigated the use of homography-based low rank models for reducing the angular dimension [START_REF] Jiang | Light field compression with homography-based low-rank approximation[END_REF], or local Gaussian mixture models in the 4D ray space [START_REF] Verhack | Steered mixture-of-experts for light field coding, depth estimation, and processing[END_REF]. The authors in [START_REF] Tabus | Lossy compression of lenslet images from plenoptic cameras combining sparse predictive coding and jpeg 2000[END_REF], use a depth based segmentation of the light field into 4D spatio-angular blocks with prediction followed by JPEG-2000.

In this paper, we consider instead geometry adapted local graph transform with a support defined by super-rays for light field compression. Super-rays have been first introduced in [START_REF] Hog | Super-rays for efficient light field processing[END_REF] as an extension of super-pixels in the 3D domain to group light rays coming from the same 3D object, i.e. to group pixels having similar color values and being close spatially in the 3D space. The method performs a k-means clustering of all light rays based on color and distance in the 3D space. Despite the limited support, computing the basis functions of the transform (i.e. the eigenvectors of the Laplacian of the graph) with a non separable spatio-angular graph (i.e. which connects pixels both within super-pixels and the ray directions), remains intractable. We therefore consider separable transforms using the optimization procedure introduced in [START_REF] Rizkallah | Geometry-aware graph transforms for light field compact representation[END_REF] to preserve angular correlation in the spatial transform domain. The local support however does not allow us to capture long term spatial dependencies, unlike efficient predictive schemes used in state of the art codecs as HEVC. To cope with this limitation, we propose a graph based spatio-angular prediction leading to a novel coding scheme for light fields.

The approach proposed in this paper is based on the observation that, when using a spatio-angular non separable transform, most of the light field energy is compacted in the low angular frequency coefficients. This motivates the design of approaches that would allow the best compression of these coefficients, e.g. by exploiting spatial correlation beyond the limits of the local graph transform support. More precisely, in this paper, we derive the equations allowing us to recover the low angular frequency coefficients of the spatio-angular transform from an encoded reference view and the high angular frequency coefficients. Thanks to these equations, the encoding of the reference view (using any classical encoder having efficient spatial predictors to exploit spatial correlation) is a way to efficiently encode the low angular frequency coefficients containing most of the light field energy. The proposed method can be seen as a graph-based prediction deriving low frequency spatio-angular coefficients from one single compressed reference view (e.g. top-left view in the experiments) and from the high angular frequency coefficients. The proposed method keeps the advantages of reduced basis function computational complexity due to the limited support, and of the energy compaction of the graph transform.

The method has been assessed in the context of quasi-lossless encoding of light fields. The compression algorithm proceeds as follows. A reference view is first compressed using a high quality encoding (we used HEVC with QP = 0 in the experiments). All the coefficients (except for the low angular frequency coefficients) of the spatio-angular graph transforms are encoded using a simple arithmetic encoder. Experimental results show that, when coupled with a powerful intra-prediction tool, the graph-based spatio-angular prediction brings a substantial gain in bitrate reaching almost 30%.

The paper is organized as follows. In section 2, a brief overview of the concept of separable graph transforms based on super-rays is provided. Section 3 gives more insight on our main contribution, the so-called graph-based spatio-angular prediction before introducing in section 4 the quasi-lossless compression scheme proposed to evaluate our prediction. Experimental results are shown in Section 5.

Separable Graph Transforms

Local Support

The first step to construct the graph transform is to define its support. Here, we consider local supports based on super-rays. To be well suited for local separable graph transforms, the super-rays need to be consistent across the views. We hence consider the approach described in [START_REF] Rizkallah | Geometry-aware graph transforms for light field compact representation[END_REF] to first compute super-pixels in the top-left view, and then project, using the top left view disparity, its segmentation labels to all the other views. The occluded pixels are assigned the label of the neighboring super-ray corresponding to the foreground objects (i.e. having the higher disparity). As for appearing pixels, they will be clustered with the background super-rays (i.e. having the lower disparity).

Spatial Graph Transforms

We define the luminance signal in a view v of the light field as x v . A graph is then defined in each super-ray k of each view v. More precisely, pixels x k,v with the same label k and belonging to the same view v are gathered in a graph

G k,v = {V k,v , E k,v },
whose edges are defined between horizontal and vertical adjacent pixels. For each graph G k,v , an adjacency matrix A k,v is defined, depicting the connections between pixels. The luminance of the pixels is considered as a real signal x k,v defined on the vertices V k,v of the graph. As proposed in [START_REF] David I Shuman | The emerging field of signal processing on graphs: Extending highdimensional data analysis to networks and other irregular domains[END_REF], a graph transform is defined based on the Laplacian matrix

L k,v = D k,v -A k,v
, where D k,v is a diagonal degree matrix whose i th diagonal element is equal to the sum of the weights of all edges incident to node i. The spatial graph transform basis is given by the eigenvectors U k,v of the Laplacian L k,v . The spatial graph transform coefficients for each graph are obtained by calculating:

xk,v = U k,v x k,v . (1)
Inversely, the luminance values of the pixels belonging to the graph are retrieved from

x k,v = U k,v xk,v . (2)
The spatial graph Laplacian in the super-ray k undergoes slight changes in different views v due to occlusions. Therefore, the basis functions of each spatial graph transform U k,v are not the same in the different views v, resulting in decreased correlation of the spatial transform coefficients xk,v across views. To circumvent this issue, we have proposed in [START_REF] Rizkallah | Geometry-aware graph transforms for light field compact representation[END_REF], a basis optimization procedure to find the basis functions U k,v that approximately diagonalize their respective Laplacians L k,v while being coherent across the views, given the scene geometry. By coherence, we mean the similarity of the eigenvectors in different views along the epipolar lines of the light field. In the following, U k,v denotes the optimized basis functions. Inside a super-ray k, the spatial transform coefficients xk,v are thus correlated between the views v. In the next section, we explain how an angular transform is used to capture this correlation.

Angular Graph Transforms

The purpose of this transform is to capture the correlation existing between the spatial transformed coefficients of each super-ray along the rays direction. For a super-ray k, we denote by spatial-band b, the vector formed by all the b th spatial transform coefficients existing in the views. More precisely, xb

k = [x k,v (b)] v=1:V, s.t. b<|x k,v |
, where V is the number of views.

We then define the graphs across views, for each band b. We link the horizontal and adjacent views together with edges with weights of 1. If we have isolated nodes, we connect them to their nearest neighbor based on the euclidean distance between view indices. Since we do not have necessarily the same number of pixels for a specific super-ray in all the views, the number of spatial coefficients after the spatial transform differs along the views. Therefore, for each spatial-band b as previously defined, we build a different graph between the views where the band b exists.

For each super-ray k, and for each band b, the Adjacency

xb k = V b k xb k . ( 3
)
3 Graph-based Spatio-Angular Prediction

Due to the high level of correlation between the different views of a light field, the energy of the transformed coefficients xb k is highly compacted in the low angular frequencies. However, some correlation remains between the coefficients of different super-rays k, and is not yet captured by the graph-transforms introduced previously. We thus propose to use already existing tools that have been proven to efficiently tract these correlations: efficient image coders involving intra prediction methods. More precisely, a reference view is coded using one of these efficient coders, and is used to recover a part of the spatio-angular spectrum at the receiver side, precisely, the one with the highest energy. In other words, we propose to compact most the light field energy on a very small set of coefficients thanks to the spatio-angular transform introduced before. And then, instead of simply coding these coefficients, we recover them from the reference view and the high angular frequencies. For that aim, we introduce a prediction tool that works as follows.

Let us consider that the view 1 is coded as a reference. For a given super-ray, its spatial transform is xk,1 = U k,1 x k,1 according to notations introduced in Sec. 2.2. For a given band b, the inverse angular transform is defined as

xb k = V b k xb k (4) ⇔      xb k (1) xb k (2) . . . xb k (N b ),      =      V b k (1, 1) V b k (1, 2) • • • V b k (1, N b ) V b k (2, 1) V b k (2, 2) • • • V b k (2, N b ) . . . . . . . . . . . . V b k (N b , 1) V b k (N b , 2) • • • V b k (N b , N b )           xb k (1) xb k (2) . . . xb k (N b )      (5)
where N b denotes the number of views where the b th band of the k th super-ray is defined.

xb k (1) = 1 V b k (1, 1)    xb k (1) -V b k (1, 2) • • • V b k (1, N b )    xb k (2) . . . xb k (N b )       . (6)
Equation ( 6) is our so-called graph-based spatio-angular prediction. The highfrequency coefficients xb k (2), . . . , xb k (N b ) can be seen as prediction coefficients, transmitted to recover the exact light field at the decoder. The basis of the linear prediction is the graph-transform basis, which makes these coefficients low-energetical and thus easy to transmit. The spatial coefficients of all the views are then retrieved from the following equation

   xb k (2) . . . xb k (N b )    =    V b k (2, 1) . . . V b k (N b , 1)    xb k (1) +    V b k (2, 2) • • • V b k (2, N b ) . . . . . . . . . V b k (N b , 2) • • • V b k (N b , N b )       xb k (2) . . . xb k (N b )   
We will show in the experimental section that the approach is very efficient when the quantization noise on the reference view is limited, hence for quasi-lossless compression. If the latter is too coarsely compressed, drift and noise amplification might appear during this prediction step. This is due to the fact that, in Equation ( 6), the prediction uses the spatial transform coefficients estimated on the reference view available at the decoder side. Further study will be dedicated to addressing this problem in the case of lossy compression.

The proposed coding scheme

Overall description of the coding scheme

The proposed quasi-lossless coding scheme is shown in Figure 1. In order to explore the inter super-rays dependencies (dependencies between different super-rays), the top left view x 1 is at first compressed with HEVC as an intra frame. It is then separated into uniform regions using the SLIC algorithm ( [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF]) to segment the image into super-pixels, and its disparity map is estimated with [START_REF] Jiang | Depth estimation with occlusion handling from a sparse set of light field views[END_REF]. The disparity values are encoded using simple arithmetic coder. Using both the segmentation map and the geometrical information, we build consistent super-rays and graphs in and across all views as explained in section 2.

We apply the separable coherent spatio angular graph transforms on the coded version of the first view (quasi-lossless coding) and the original versions of all other frames to compact their energy in fewer coefficients. Since the first view is coded with very small QP, we are then almost sure that we are not adding angular incoherence between the first view and the rest of the views. Once we have the spatio-angular coefficients, instead of sending the whole spectrum with simple arithmetic coding, we propose to make use of our graph-based prediction and therefore deriving low frequency spatio-angular coefficients in the decoder side from the top-left reference view coded with HEVC-Intra and from the high angular frequency coefficients.

We thus send, for all super-rays, the AC coefficients of the N k,1 first bands after the second angular graph transform. (N k,1 is the number of pixels belonging to the super-ray k in view 1). We come back to the detailed coding of those coefficients in the sequel. Note that to construct the spatial and angular graphs on the decoder side, we do not need extra information. Since the decoder already receives the top left image and the disparity values, with the SLIC algorithm, it can deduce the segmentation map and the super-ray segmentation used to for local graph transforms design and application. Furthermore, the decoder have received the spatio-angular high frequency coefficients. Using those two kinds of information, it can predict the DC spatioangular components obtained after the angular graph transform and then reconstruct the luminance values of all the views as explained in Section 3.

The coding of the spatio-angular high frequency coefficients

After applying the spatio-angular graph transforms on the whole super-ray, all frequency coefficients are grouped into a three-dimensional array y where y(k, b, v) is the v th transformed coefficient of the band b for the super-ray k. Using the observations on all the super-rays in a training dataset (Rock ), we can find the best ordering for quantization. We first sort the variances of coefficients with enough observations in decreasing order. Using this scanning order, we assign a class number to each super-ray. More precisely, each super-ray with N coefficients belongs to class i if the mean energy per high frequency coefficient is less than 1, where the high frequency coefficients are the last round(N × i/4) coefficients following the scanning order of the super-rays coefficients defined previously.

We start by finding the super-rays in the first class than remove them from the search space before finding the other classes, and idem for the following steps. We code a flag with an arithmetic coder to give the information of the class of super-rays to the decoder side. In class i, the last round(N × i/4) coefficients of each super-ray are discarded. Then, the rest of the high frequency spatio-angular coefficients are quantized uniformly with a small step size Q = 1. Then, they are grouped in 32 uniform groups to enter a simple arithmetic coding assuming an iid model. To evaluate our proposed spatio-angular prediction method in the full coding scheme, we apply it on real light fields captured by plenoptic cameras from the datasets used in [START_REF] Khademi Kalantari | Learning-based view synthesis for light field cameras[END_REF] and [START_REF] Viola | A graph learning approach for light field image compression[END_REF]. To avoid the strong vignetting and distortion problems on the views at the periphery of the light field, we only consider the 8 × 8 central subaperture images cropped to 364×524 in [START_REF] Khademi Kalantari | Learning-based view synthesis for light field cameras[END_REF], and 9×9 cropped to 432×624 from [START_REF] Viola | A graph learning approach for light field image compression[END_REF]. Some of the light fields considered are shown in Figure 2. The full set of light fields considered for the test is: Flower2, Cars ,Rock and Seahorse from the dataset in [START_REF] Khademi Kalantari | Learning-based view synthesis for light field cameras[END_REF] and StonePillarInside and Friends from the dataset of ICIP challenge 2017 and used in [START_REF] Viola | A graph learning approach for light field image compression[END_REF]. The method used to estimate the disparity of the top-left views is described in [START_REF] Jiang | Depth estimation with occlusion handling from a sparse set of light field views[END_REF]. Examples of the disparity maps provided are shown in Figure 3. A sparse set of disparity values and the segmentation maps computed with SLIC [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF] are used to construct local graph supports as described in Section 2. In our experiments, we fix the number of super-rays for all datasets to 4000. Indeed, an ideal segmentation for compression purposes, would be the result of a local rate-distortion optimization (as previously proposed in [START_REF] Su | Rate-distortion optimized super-ray merging for light field compression[END_REF]). For the sake of simplicity in this paper, we use a simple non adaptive segmentation.

The importance of exploring the long term dependencies between super-rays

As explained before, we aim at compacting the energy of the light field luminance in fewer coefficients, then we will predict them based on the knowledge of a coded topleft view and the high frequency components that most likely are sparse. To illustrate the effect of our spatio-angular transforms on the energy compaction, we compute in Table 1 the percentage of total energy that resides in the predicted DC spatio-angular bands xb k (1) ∀k, b. Those bands are actually not sent but predicted on the decoder side, and only the energy that is left is coded as extra information to reconstruct the original light field. We can observe that most of the energy is compacted in those coefficients, which shows the efficiency of the spatio-angular decorrelation of the Graph transforms. Thanks to the equations introduced in Section 3, the encoding of the top-left view (using any classical encoder having efficient spatial predictors to exploit spatial correlations) can be seen as a way to efficiently encode those DC spatio-angular frequency coefficients which contain most of the light field energy. We choose the HEVC intra in order to exploit long term dependencies between super-rays. More precisely, we code the top left view with HEVC Intra. Then in Table 2, we compare its performance with respect to using a simple arithmetic coder for the DC spatioangular frequency coefficients. The QP of HEVC is set to 0. As for the arithmetic coder, the observations of each DC spatio-angular frequency coefficient xb k (1) ∀k, b are coded separately. We can see for most of the datasets, the gain obtained with HEVC intra when compared to arithmetic coding without any spatial prediction.

Comparative assessment

To assess our Graph-based Spatio-Angular Prediction, we evaluate it in the context of the quasi-lossless coding scheme in Section 4 against a complete HEVC based scheme with a QP set to 0 and a GOP of 4. The HEVC version used in the tests is HM-16.10. The light fields are coded following a raster scanning starting with the top-left view as a reference intra-coded frame. Results are reported in Table 3 where we compare mainly the rate needed to code a Light Field under quasi-lossless settings (based on visual quality assessment, we consider a PSNR higher that 50 dB as a quasi-lossless compression). A substantial gain in bitrate is observed while preserving a high quality of the reconstructed light fields. This can be justified by the efficiency of our spatio-angular graph transforms in of compaction along with the ability of HEVC-intra to effectively exploit spatial correlations in the reference view.

Conclusion

In this paper, we have explored local separable spatio-angular graph transforms for light fields compact representation. The limited support of local transforms may not allow us to exploit long term spatial dependencies. To cope with this limitation, we have proposed a novel approach to leverage the good spatial de-correlation properties of traditional codecs (e.g. HEVC intra), making use of efficient predictors, into local spatio-angular graph transforms. A reference view coded with any efficient codec is used to predict low angular frequency transform coefficients that, together with the transmitted high angular transform coefficients, allow recovering the entire graphbased representation. The scheme has been assessed for high quality (quasi-lossless) coding of light fields with narrow baselines. Note that the prediction mechanism, in case of coarser quantization, tends to amplify the quantization noise on the reconstructed data representation and the method need to be tuned for wider baselines. Those issues are left for further study.

  A b k and degree D b k matrices are used to compute the inter-view Laplacian as L b k = D b k -A b k . The eigenvectors V b k of L b k are then used to define the angular graph transform basis. The angular transform coefficients are obtained by calculating:

Figure 1 :

 1 Figure 1: Overview of proposed coding scheme.

Figure 2 :

 2 Figure 2: An example set of light fields used in our experiments. Only the top-left view is shown for illustration purpose. From left to right: Flower2, Cars, StonePillarsInside and Friends.

Figure 3 :

 3 Figure 3: An example set of top-left view disparity maps used in our experiments

Table 1 :

 1 Percentage of energy residing in the DC spatio-angular bands xb k (1) ∀ k, b

	Light Fields	Energy Percentage in xb k (1) ∀ k, b
	Flower 2	99.02 %
	Cars	99.34 %
	Rock	98.45 %
	Seahorse	98.73 %
	Stone Pillars Inside	98.26 %
	Friends	99.80 %

Table 2 :

 2 Rate comparison between HEVC intra-coding of the first view and entropy coding of all the DC spatio-angular bands xb

	Light Fields	HEVC intra Entropy coding
	Flower 2	1.0 Mbits	1.48 Mbits
	Cars	1.06 Mbits	1.54 Mbits
	Rock	1.02 Mbits	1.38 Mbits
	Seahorse	0.72 Mbits	1.22 Mbits
	Stone Pillars Inside 1.81 Mbits	1.66 Mbits
	Friends	1.62 Mbits	2.04 Mbits

k (1) ∀ k, b using an arithmetic coder assuming an iid model

Table 3 :

 3 Rate comparison between our proposed approach and a traditional HEVC-inter according to a raster scan order for high PSNR (> 50 dB) (54.2033 dB) 2.4799 bpp (55.1969 dB) Cars 3.6688 bpp (54.0812 dB) 2.6258 bpp (55.2009 dB) Rock 3.2700 bpp (53.7601 dB) 2.0162 bpp (54.7765 dB) Seahorse 2.4751 bpp (54.3804 db) 1.9762 bpp (55.2844 dB) Stone Pillars Inside 4.9017 bpp (52.1036 dB) 3.3094 bpp (55.0022 dB) Friends 3.5400 bpp (52.7986 dB) 2.4436 bpp (54.8196 dB)

	Light Fields	HEVC Raster Scan	Ours
	Flower 2	3.3129 bpp

This work has been supported by the EU H2020 Research and Innovation Programme under grant agreement No 694122 (ERC advanced grant CLIM).


Pierre Allain

Laurent Guillo

Christine Guillemot

LIGHT FIELD DENOISING USING 4D ANISOTROPIC DIFFUSION

Keywords: Light Field, regularization, anisotropic diffusion, denoising

In this paper, we present a novel light field denoising algorithm using a vector-valued regularization operating in the 4D ray space. More precisely, the method performs a PDE-based anisotropic diffusion along directions defined by local structures in the 4D ray space. It does not require prior estimation of disparity maps. The local structures in the 4D light field are extracted using a 4D tensor structure. The paper then describes the strategy retained for setting the diffusion tensor parameters for the targeted denoising application. It then analyzes the influence of the model parameters on the denoising performance. Experimental results show that the proposed denoising algorithm performs well compared to state of the art methods while keeping tractable complexity.

INTRODUCTION

Light field imaging has emerged as a promising technology for a variety of applications including photorealistic rendering, computational photography and computer vision. Various camera designs have been proposed to capture real light fields going from camera arrays [START_REF] Wilburn | High performance imaging using large camera arrays[END_REF], to single cameras mounted on moving gantries, and plenoptic cameras [START_REF] Ng | Light Field Photography[END_REF], [START_REF] Georgiev | Superresolution with the focused plenoptic camera[END_REF]. These camera architectures present different characteristics in terms of spatial and angular resolution. Plenoptic cameras rely on multiplexing and use an array of microlenses placed in front of the sensor to capture multiple low resolution views in one 2D sensor image. This easy way to capture multiple view-points, suffers from a low spatial resolution and reduced light throughput, resulting in visible noise artefacts, compared to 2D classical cameras.

This paper addresses the problem of denoising in dense light fields. Denoising is an ill-posed problem which, to be solved, requires introducing some prior knowledge on the kind of typical images we try to restore, which helps restricting the class of admissible solutions. The authors in [START_REF] Mitra | Light field denoising, light field superresolution and stereo camera based refocussing using a gmm light field patch prior[END_REF] use a Gaussian mixture model (GMM) prior for light field patches, and then exploit this prior within a Bayesian inference framework for light field denoising. In [START_REF] Moghaddam | Light field denoising: exploiting the redundancy of an epipolar sequence representation[END_REF], the light field is denoised by stacking EPIs in a 3D volume using the VBM4D video denoising method [START_REF] Boracchi | Video denoising, deblocking, and enhancement through separable 4-d nonlocal spatiotem-poral transforms[END_REF], while the authors in [START_REF] Dansereau | Light field image denoising using a linear 4d frequency-hyperfan all-in-focus filter[END_REF] apply a hyper-fan shaped filter in a 4D discrete Fourier transform domain. The BM3D denoising algorithm [START_REF] Dabov | Image denoising by sparse 3-d transform-domain collaborative filtering[END_REF] is extended in [START_REF] Alain | Light field denoising by sparse 5d transform domain collaborative filtering[END_REF] to light fields.

. This work was supported by the EU H2020 Research and Innovation Programme under grant agreement No 694122 (ERC advanced grant CLIM) Smoothness priors have been widely used together with regularization schemes for inverse problems in 2D imaging such as denoising. The only approach of regularization proposed so far for light field denoising is the one of [START_REF] Goldluecke | The variational structure of disparity and regularization of 4d light fields[END_REF] using a global variational framework integrating contributions of regularizers on vector fields computed in EPI 2D spaces. The authors use 2D structure tensors computed in epipolar plane images (EPI) to estimate disparity between light field views. With the help of the estimated disparity, they project all the views on a target viewpoint and solve the regularization using a global optimization framework with total variation priors. The final regularizer is written as the contributions of regularizers on vector fields in EPI 2D spaces, with user-defined constants which adjust the amount of smoothing on the different views or EPIs. Errors in disparity maps may however translate into annoying artefacts in particular along edges.

In this paper, we propose instead a regularization method operating directly in the 4D ray space that does not require prior estimation of disparity maps. The method performs a PDE-based diffusion with anisotropy steered by a tensor field based on local structures in the 4D ray space that we extract using a light field 4D tensor structure. To enhance coherent structures, the smoothing along directions, surfaces, or volumes in the 4D ray space is done with a diffusivity coefficient that is derived from the amount of local variations in the 4D space. More precisely, the diffusivity coefficient is computed as a function of the 4 eigenvalues of the 4D structure tensor. The diffusion is isotropic in homogenous regions with no strong structures, whereas it follows structures in non homogenous regions.

The proposed regularization model has been applied for denoising light fields from the EPFL dataset [START_REF] Rerabek | New light field image dataset[END_REF]. Experimental results show that the proposed method compares well to the best state of the art method [START_REF] Alain | Light field denoising by sparse 5d transform domain collaborative filtering[END_REF], even outperforming it at high noise levels, while being considerably faster. We refer to our method as Anisotropic Diffusion Denoising 4D (ADD4D).

4D ANISOTROPIC DIFFUSION

Let L(x, y, u, v) denote the 4D representation of a light field, describing the radiance of a light ray parameterized by its intersection with two parallel planes [START_REF] Gortler | The lumigraph[END_REF], and where (u, v) denote the angular (view) coordinates and (x, y) the spatial (pixel) coordinates.

4D structure tensor

To perform light field regularization along pre-defined directions which would preserve edges, one has to first characterize local structures in the image, or here the light field. While light field regularization methods are usually applied per Epipolar Plane Image (EPI) of the light field with directions given by 2D structure tensors often smoothed with a 2D Gaussian kernel, we consider instead a direct computation of 4D structure tensors in the 4D ray space, which is defined as the symmetric and semipositive definite 4 × 4 matrix [13]:

T = 1 C C c=1 ∇L c ⊗ ∇L c , ( 1
)
where C is the number of color channels. The notation ∇L c denotes the gradient or partial derivatives in the 4 dimensions of the color channel L c . A smoothed version T ω = T * G ω (x, y, u, v) of the structure tensor is in practice computed to retrieve a more coherent geometry, where G ω (x, y, u, v) is a Gaussian kernel smoothing along the four light field dimensions. Its spectral elements, i.e. the eigenvectors ν i , i = 1 . . . 4 and the eigenvalues λ i , i = 1 . . . 4 of the 4D structure tensor T respectively give the orientations of the local structures in the 4D rays space and the vector-valued variations along these structures. We therefore have:

T ω = 4 i=1 λ i ν i ν T i (2)
The eigenvector ν 1 with the largest eigenvalue λ 1 gives the dominant orientation of the local structures in the 4D ray space. The ratios between the eigenvalues gives information on the type of structure present in the local neighborhood.

TRACE-BASED PDE DIFFUSION

We consider a trace-based PDE as introduced in [START_REF] Perona | Scale space and edge detection using anisotropic diffusion[END_REF] for 2D images, which we here extend to 4D light fields as:

∂L c ∂t = trace(DH c ) (3)
where H c is the 4D Hessian matrix of the c channel of L, defined as:

H c =       ∂ 2 Lc ∂x 2 ∂ 2 Lc ∂x∂y ∂ 2 Lc ∂x∂u ∂ 2 Lc ∂x∂v ∂ 2 Lc ∂y∂x ∂ 2 Lc ∂y 2 ∂ 2 Lc ∂y∂u ∂ 2 Lc ∂y∂v ∂ 2 Lc ∂u∂x ∂ 2 Lc ∂u∂y ∂ 2 Lc ∂u 2 ∂ 2 Lc ∂u∂v ∂ 2 Lc ∂v∂x ∂ 2 Lc ∂v∂y ∂ 2 Lc ∂v∂u ∂ 2 Lc ∂v 2       (4)
and D is a diffusion tensor to be defined according to the regularization goals.

To be adapted to the local structures in the 4D ray space of the light field, the diffusion process must be driven by a specific tensor, called diffusion tensor chosen so that it has the same set of eigenvectors ν i , i = 1 . . . 4 as the structure tensor T. A diffusion tensor is therefore constructed as:

D = 4 i=1 α i ν i ν T i , (5)
where the coefficient α i set the smoothing strength along the direction given by ν i , and are chosen according to the goal of the diffusion.

Setting the diffusion tensor parameters

For denoising purpose, the coefficients α i controlling the strength of the diffusion have to depend on the quantity of local variation. Considering ψ = 4 i=1 λ i as a measure of local variation [START_REF] Tschumperlé | Anisotropic Diffusion Partial Differential Equations in Multi-Channel Image Processing : Framework and Applications[END_REF], we set these parameters as α i = f i (ψ). The functions f i (ψ) are chosen to abide by denoising goal. Firstly, the diffusion should be isotropic in regions where there is no apparent direction, i.e. where L is homogeneous:

lim D ψ→0 = aI ⇒ f 1 (0) = f 2 (0) = f 3 (0) = f 4 (0). (6)
Conversely, in regions where structures are clearly appearing, no diffusion must be performed in order to preserve geometry, that is:

lim D ψ→inf = 0 ⇒ f 1 (ψ) = f 2 (ψ) = f 3 (ψ) = f 4 (ψ) = 0. (7)
Also, the diffusion must be performed along the orientations of the local structures, given by the eigenvectors having the lowest eigenvalues. We assume that more structure there is, more the eigenvectors corresponding to the lowest eigenvalues must drive the diffusion. For instance, for a high quantity of structure, the diffusion should only be performed in the ν 4 direction. In other words, if the structure is clearly shaped, the diffusion should only correct structure singularities with the highest degree. Applying such reasoning to the other eigenvectors, we obtain the condition:

f 1 (ψ) ≤ f 2 (ψ) ≤ f 3 (ψ) ≤ f 4 (ψ). (8)
We therefore propose the function f i meeting the conditions of equations 6, 7 and 8:

f i (ψ) = e -ψ ρ i γ , with ρ 1 ≤ ρ 2 ≤ ρ 3 ≤ ρ 4 (9)
where ρ i is a parameter related to the quantity ψ of local variations corresponding to a point where the diffusion strength f (ψ) significantly varies. The parameter γ defines the exponential behaviour of the function. In the particular case where γ = 2, the parameter ρ i corresponds to the standard deviation of the function f (ψ). This parameterization allows the tuning of the diffusion process on the basis of structure properties.

Model tuning

We discuss here the influence of the model parameters for denoising. Instead of parameterizing independently every ρ i , we propose the following expression:

ρ i = ρ 0 α i-1 , ( 10
)
which allows a simplification of the model to two parameters. Figure 3 plots different functions f i (ψ) parameterized by ρ i . We set α and γ to 2 and study the influence of ρ 0 on the denoising process. One can note that α is related to an anisotropic ratio, while ρ 0 is related to a local variation intensity. We present denoising results with different set-ups for a single light field. The model integration is performed using the same time step of 0.1 to raise parameters influence of the denoising. Fig. and2b, show that, for low and medium noise (σ = 10 255 , σ = 40 255 ), a low value of ρ 0 (respectively 0.01 and 0.03) prevents the model to damage the light field structures. On the other hand, a low value of ρ 0 is not suitable for high noise levels as shown in 2d, where the model is unable to diffuse along the structures, hence many iterations are needed to remove the noise. For high values of ρ 0 , we can observe a more efficient behaviour regarding noise reduction. But in this case, the model can damage the structures (more than with low values of ρ 0 ), as we can see when the number of iterations increases.

Fig. 2d-f show that, for the different values of the parameter ρ 0 , when increasing the number of iterations, the model converges to the same PSNR curve, at both low and high noise level. This shows that, even for very different intensities of added noise, the model is mostly responsible for the new structures of the light field. In other words, noise is simply cleaned more or less quickly while the model is on its way to its intrinsic regularization evolution. [START_REF] Alain | Light field denoising by sparse 5d transform domain collaborative filtering[END_REF]). PSNR in dB.

DENOISING RESULTS

The proposed denoising model is applied on 12 light fields of the EPFL dataset [START_REF] Rerabek | New light field image dataset[END_REF] extracted using Dansereau's Matlab light field toolbox [START_REF] Dansereau | Decoding, calibration and rectification for lenselet-based plenoptic cameras[END_REF] in the same conditions as in [START_REF] Alain | Light field denoising by sparse 5d transform domain collaborative filtering[END_REF]. They are chosen to allow comparison to state of the art methods. Because of the strong irregularities in the angular dimension resulting from the vignetting effect, we chose to denoise only the subapertures that are not affected by the later and therefore limit the size of light fields to 11×11. Indeed, because of its regularization properties, our method would also correct artefacts at angular corners (i.e. the vignetting effect), hence the measured PSNR would not only reflect the denoising performance.

For our experiments we chose to use the same model parameterized by ρ 0 = 0.05, α = γ = 2, and tensor smoothing is set to ω = 5 in the spatial and angular directions. Integration of Eq. 3 is performed using the Runge-Kutta 4 scheme, and the time step is set to 0.1. One iteration takes 7 minutes to compute on a i7-6600U without any parallelization, and best PSNR is usually obtained from 4 iterations to 20 depending on the light field and the model. For comparison, LFBM5D performs approximately in 7 hours on an octo core processor (but without angular cropping). On a 10 iterations basis for our method, this gives an estimated 26 times faster result. Detailed results are provided 1 to show angular consistency. In figure 4, we observe that our model is able to remove low noise while keeping fine original structures. Our model is also able to deal with high level of noise as shown in figures 1 and 5, where we can see in the last that the result qualitatively matches the best existing method. In table 1 we can see the denoising model we used outperforms other methods for a high level of noise, while comparing well for low levels.

CONCLUSION

In this paper we proposed a new approach for denoising light fields. Using structural properties in the four dimension space, we showed that tensor oriented diffusion can drastically reduce noise while keeping most of the structures of the subaperture images. We believe that this work can be extended with different denoising models depending on scene's disparities and structures, and can highly benefit of parallelization.

Fig. 1 :

 1 Fig. 1: 4D tensor driven diffusion is able to denoise highly damaged light fields. Center views for Vespa and DangerDeMort from the EPFL dataset. (a,c) With Additive Gaussian noise (σ = 100 255 , PSNR= 8.13 dB); (b,d) After denoising giving global PSNR on all subapertures of 29.31 dB (b) and 28.50 dB (d).

Fig. 3 :

 3 Fig.3: Functions f i (ψ) controlling smoothing strength of the diffusion tensor along orientations defined by the four eigenvectors ν i , with ρ 0 = 0.05.

10 Fig. 2 :Fig. 4 :Fig. 5 :

 10245 Fig. 2: Denoising behaviour for different noise levels applied to the Bike light field. Each curve represents the PSNR evolution along the iterations of the diffusion for: (top row) -different values of the parameter ρ 0 of the model, and for different noise levels For low and medium noise (a,b), we can see that the best PSNR is obtained for a low value of the parameter σ 0 ; (bottom row) -different noise levels, and for different values of the parameter ρ 0 of the model .

Table 1 :

 1 Comparison of light field denoising methods (state of the art results are courtesy of

	Method	σ = 10 255 σ = 20 255 σ = 30 255 σ = 40 255 σ = 50 255
	HF4D [7]	31.07	25.80	22.61	20.34	18.59
	BM3D [8]	35.42	32.85	31.36	30.25	29.32
	BM3D Epi [17]	36.09	33.48	31.90	30.71	29.67
	VBM4D [18]	36.07	33.52	31.92	30.67	29.63
	VBM4D Epi [5]	36.13	33.51	31.92	30.72	29.72
	LFMB5D [9]	36.50	34.21	32.87	31.84	30.99
	ADD4D	36.13	33.88	32.67	31.84	31.19

, PSNR 14.15 dB. (b) Result obtained using LFBM5D, PSNR 29.91 dB. (c) Using our method, 19 iterations, PSNR 30.10 dB.


Xiaoran Jiang

Jinglei Shi

Christine Guillemot

A LEARNING BASED DEPTH ESTIMATION FRAMEWORK FOR 4D DENSELY AND SPARSELY SAMPLED LIGHT FIELDS

Keywords: depth estimation, light field, deep learning, multi-view stereo

This paper proposes a learning based solution to disparity (depth) estimation for either densely or sparsely sampled light fields. Disparity between stereo pairs among a sparse subset of anchor views is first estimated by a fine-tuned FlowNet 2.0 network adapted to disparity prediction task. These coarse estimates are fused by exploiting the photo-consistency warping error, and refined by a Multi-view Stereo Refinement Network (MSRNet). The propagation of disparity from anchor viewpoints towards other viewpoints is performed by an occlusion-aware soft 3D reconstruction method. The experiments show that, both for dense and sparse light fields, our algorithm outperforms significantly the state-of-the-art algorithms, especially for subpixel accuracy.

INTRODUCTION

Unlike traditional cameras, light field (LF) devices capture ray intensities emitted by a 3D point along different orientations. Due to this rich description of the 3D scenes, great attention has been given to efficient and robust algorithms for depth estimation exploiting LFs. According to different kinds of images that the estimation depends on, existing algorithms can be classified into several main categories: methods based on Sub-Aperture Images (SAI) [START_REF] Heber | Shape from light field meets robust PCA[END_REF][START_REF] Jeon | Accurate depth map estimation from a lenslet light field camera[END_REF][START_REF] Chen | Fast and accurate optical flow based depth map estimation from light fields[END_REF], on Epipolar Plane Images (EPI) [START_REF] Wanner | Variational light field analysis for disparity estimation and super-resolution[END_REF][START_REF] Zhang | Robust depth estimation for light field via spinning parallelogram operator[END_REF][START_REF] Johannsen | What sparse light field coding reveals about scene structure[END_REF][START_REF] Shin | Epinet: A fully-convolutional neural network using epipolar geometry for depth from light field images[END_REF] or on refocused images [START_REF] Tao | Depth from combining defocus and correspondence using light-field cameras[END_REF][START_REF] Ting-Chun | Occlusion-aware depth estimation using light-field cameras[END_REF]. However, most of these methods are designed for dense view sampling. Although some algorithms only use a subset of views, e.g. horizontal or crosshair viewpoints, disparity between the exploited views should remain small.

On the contrary, the authors in [START_REF] Huang | Empirical bayesian light-field stereo matching by robust pseudo random field modeling[END_REF] employ an empirical Bayesian framework, which adapts parameters according to scene statistics. This algorithm is free of additional cues exploiting dense view sampling, e.g. phase shift [START_REF] Jeon | Accurate depth map estimation from a lenslet light field camera[END_REF], spinning parallelogram operator [START_REF] Zhang | Robust depth estimation for light field via spinning parallelogram operator[END_REF], defocus cue [START_REF] Tao | Depth from combining defocus and correspondence using light-field cameras[END_REF] and structured tensor [START_REF] Wanner | Variational light field analysis for disparity estimation and super-resolution[END_REF], thus it is robust and relevant for both dense and sparse LFs. In [START_REF] Jiang | Depth estimation with occlusion handling from a sparse set of light field views[END_REF], densely sampled depth (disparity) information is inferred from a subset of sparsely sampled light field views. Disparity between sparse views are first estimated and refined using optical flow estimator and edge-preserving filtering, then they are propagated to other viewpoints by exploiting angular correlation.

This work has been funded by the EU H2020 Research and Innovation Programme under grant agreement No 694122 (ERC advanced grant CLIM).

Meanwhile, deep learning has met great success in binocular vision, and similarly in optical flow estimation. As one of the pioneers, FlowNet [START_REF] Dosovitskiy | Flownet: Learning optical flow with convolutional networks[END_REF] employs an end-to-end encoder-decoder architecture with additional skip-connections between contracting and expanding parts. Being a variant of FlowNet, DispNet [START_REF] Mayer | A large dataset to train convolutional networks for disparity, optical flow, and scene flow estimation[END_REF] considers 1D correlation instead of 2D correlation to better adapt to disparity estimation task. By stacking several elementary networks, each of them being similar to FlowNet, FlowNet 2.0 [START_REF] Ilg | Flownet 2.0: Evolution of optical flow estimation with deep networks[END_REF] significantly improves the prediction accuracy. A cascade framework is proposed in [START_REF] Pang | Cascade residual learning: A two-stage convolutional neural network for stereo matching[END_REF], which corrects the disparity initialization by learning in a supervised fashion residual signals across multiple scales. Deep learning has been also successfully applied to light field depth estimation. Among them, EpiNet [START_REF] Shin | Epinet: A fully-convolutional neural network using epipolar geometry for depth from light field images[END_REF] achieves the state-of-the-art performance by using a multi-stream network, each stream exploiting one angular direction of light field views: horizontal, vertical, left or right diagonal directions. But this approach is well suited for dense light fields only.

In this paper, we focus on how to handle either a dense or a sparse light field for disparity/depth estimation. Similar to [START_REF] Jiang | Depth estimation with occlusion handling from a sparse set of light field views[END_REF], our algorithm only exploits a sparse subset of anchor views (four corner views), and generates one disparity map for every viewpoint of the light field. Multi-view stereo (MVS) is implemented in a deep learning based cascaded framework. A pre-trained FlowNet 2.0 is fine-tuned by pairs of stereo images, and the obtained model is used to estimate disparity between pairs of anchor views, arranged horizontally or vertically. These coarse estimates are then fused at each anchor viewpoint by exploiting the warping error from other anchor viewpoints, and then refined by a second convolutional neural network (CNN), which we call Multi-view Stereo Refinement Network (MSRNet). For better subpixel accuracy of the disparity values, views are upsampled before being fed to CNNs. Correspondingly, the output disparity maps are rescaled. The propagation of disparity from anchor viewpoints towards other viewpoints is performed by an occlusion-aware soft 3D reconstruction method.

LEARNING-BASED DISPARITY ESTIMATION

2.1. Finetuned FlowNet 2.0 for stereo We take the four corner views I 1,1 , I 1,V , I U,1 and I U,V (c.f. Fig. 1) as the set of anchor light field views. Arguably, these distinct views on the extreme corners of a densely sampled light field contain all color and geometric information, from which the whole light field can be reconstructed [START_REF] Jiang | Light field compression using depth image based view synthesis[END_REF]. Disparity is estimated between pairs of anchor views. To work within the stereo vision framework, we assume that the light field is well rectified, i.e., a scene point moves only horizontally from position (1, 1) to (U, 1), and only vertically from position [START_REF] Heber | Shape from light field meets robust PCA[END_REF][START_REF] Heber | Shape from light field meets robust PCA[END_REF] to (1, V ). In order to apply a single network to deal with both horizontal and vertical image pairs, vertical pairs are rotated 90 degrees. During the training process, their corresponding ground truth disparity maps are also rotated 90 degrees. For each anchor view Ir, two coarse estimates d r and d r are obtained, respectively using horizontal and vertical pairs.

Multi-view disparity refinement

In low-level vision tasks, such as segmentation, denoising and optical flow estimation, results can be generally improved with post-facto refinement. Similar to [START_REF] Pang | Cascade residual learning: A two-stage convolutional neural network for stereo matching[END_REF], our Multi-view Stereo Refinement Network (MSRNet) takes as input the coarse disparity estimates from a first stage network (finetuned FlowNet 2.0), and the refinement of these maps is implemented by a multi-scale encoder-decoder structure. Differently, our scheme generalizes the binocular vision refinement task to MVS scenario, and our network is flexible for any number of initialized disparity maps.

Disparity fusion

In order to fuse multiple coarse disparity estimates to a single one, we propose to leverage photo-consistency warping errors.

I i r (d r )
denotes the warped image from position i to r (both i and r are corner positions) by using disparity d r . Traditionally, warping errors from different corner viewpoints are aggregated by computing the sum or average of them:

Enocc(d r ) = mean i E(Ir, I i r (d r )) , (1)
with E(I, I ) the pixel-wise sum of square errors for the three color components between image I and I :

E(I, I ) = C∈{R,G,B} (I C -I C ) 2 .
(2)

The term Enocc measures accurately the pixel warping error in occlusion-free areas, and therefore reflects well the disparity accuracy at these pixels, but fails in occlusion zones. In fact, a high error in occlusion zones is due to interpolation in large holes, rather than to disparity inaccuracy.

Compared to stereo vision, multi-view stereo should provide a better modeling for occlusion, since a pixel occluded from one viewpoint may be viewed from another viewpoint. Therefore, we define a second error map Eocc: (

) 3
In general, error due to interpolation is smaller for non-occluded pixels than occluded ones. At a pixel p that can be seen in the warped view I i r , but not in I i r for any other corner viewpoint i, the error Eocc(d r , p) may equals to the value E(Ir, I i r (d r )) at pixel p, which more faithfully reflects the disparity accuracy for this pixel.

Hence, disparity inaccuracy is better modeled by Enocc in occlusion-free areas, and by Eocc in non-overlapped occlusion zones. Two fused disparity maps can be obtained: at each pixel p, one disparity value is selected among the candidate values (5)

One may also approximate the binary occlusion mask M by identifying the pixels that correspond to high values of Enocc:

M (p) = 1, if min Enocc(d r , p), Enocc(d r , p) > θ 0, otherwise, (6)
the threshold θ being chosen as the top 90% Enocc value. Compared to [START_REF] Jiang | Depth estimation with occlusion handling from a sparse set of light field views[END_REF] where the occlusion mask is explicitly computed based on disparity values, Eq. ( 6) shows to be a good approximation and above all more time efficient, which is compatible to the learning-based framework. Finally, the unique disparity value per pixel at corner viewpoint r is computed as:

Dr(p) = D nocc r (p) × 1 -M (p) + D occ r (p) × M (p). (7)

MSRNet

A cascade residual learning framework is proposed in [START_REF] Pang | Cascade residual learning: A two-stage convolutional neural network for stereo matching[END_REF] to refine disparity for binocular vision: residuals of disparity are learned explicitly by an encoder-decoder multi-scale network, which is able to correct the disparity initialized by a first stage CNN. Residual signals are supervised at each resolution scale.

The input of this refinement network is a set of 5 images: the left image I L , the right image I R , the initialized disparity map, warped image ĨL and warping error. In a MVS scenario, it is obvious that this scheme is no longer applicable: the multiplication of the number of stereo pairs, as well as the number of initialized disparity maps, will rapidly enlarge the size of the network and make the learning inefficient. Moreover, this scheme cannot easily adapt to the varying number of input images. This is why we chose to fuse the horizontal and vertical disparity estimates (c.f. Section 2.2.1) before residual learning. Regardless of potentially different numbers of stereo pairs, only one fused disparity map Dr will be obtained at corner position r. Two other images are also fed to MSRNet: the color image Ir as guidance, and a binary map B indicating unreliable pixels. These pixels are supposed to be located near object contours which are determined using canny edge detectors. We detect the contours on the disparity map Dr that do not correspond to those on the color image Ir. The binary mask is therefore computed as The rest of the network structure is similar to the one of [START_REF] Pang | Cascade residual learning: A two-stage convolutional neural network for stereo matching[END_REF]. Compared to [START_REF] Pang | Cascade residual learning: A two-stage convolutional neural network for stereo matching[END_REF], we add a gradient term L 2 in the loss function, which leads to smoother maps:

L = λ 1 L 1 (D, D gt ) + λ 2 L 2 (D, D gt ) (8)
where L 1 is sum of absolute differences (SAD), and the term L 2 is defined as the sum of

DISPARITY PROPAGATION

The obtained disparity maps D1,1 , D1,V , DU,1 and DU,V at the 4 corners should be propagated to other viewpoints where we assume that the color information is absent. In [START_REF] Jiang | Depth estimation with occlusion handling from a sparse set of light field views[END_REF], disparity values are projected to novel view positions, and low rank assumption of the matrix containing all warped maps is exploited to inpaint missing disparity values for occluded pixels. However, for light fields with wide baseline, this low rank assumption cannot be held anymore.

Authors in [START_REF] Penner | Soft 3d reconstruction for view synthesis[END_REF] have proposed a soft 3D reconstruction method for light field view synthesis. We apply this idea to disparity propagation. For each corner view, a consensus volume is constructed, each voxel encoding the consensus score of all voters (4 corner view disparity estimates) for the existence of a surface with a certain discretized disparity/depth value at a pixel p. To compute the consensus volume, a vote volume and a confidence volume are previously constructed, half of the voters with low confidence values abstain to vote. Occlusion/visibility is modeled by the soft visibility function: a pixel can be viewed at a depth z only if the consensus for a surface existing at any depth z < z is low (note that the disparity is inversely proportional to the depth). These soft visibility values at corner viewpoints are then used as weights for merging warped disparity maps in other viewpoints.

EXPERIMENTS

Datasets

In order to train our deep neural networks, two synthetic LF datasets with different disparity ranges have been created using Blender software [START_REF]Blender software[END_REF]. 3D models, available under CC0 or CC license, have been downloaded from Chocofur [START_REF]Chocofur 3d models website[END_REF] and Sketchfab websites [START_REF]Sketchfab 3d models website[END_REF]. The sparse light field (SLF) dataset contains scenes with disparity range [-20, 20] between adjacent views, whereas the disparity range is [-4, 4] for the dense light field (DLF) dataset. Each LF has the same resolution 512×512×9×9. We provide for each view of the LFs, a color view, a ground truth depth map and its corresponding disparity map.

Among the 53 scenes in the SLF dataset, 44 scenes are served as training data, and 9 others as valid data. For the DLF dataset, the training set contains 38 scenes and the valid set contains 5 scenes. To our knowledge, they are among the first synthetic light field datasets providing depth information for every viewpoint of the light fields (HCI 4D light field dataset [START_REF] Honauer | A dataset and evaluation methodology for depth estimation on 4d light fields[END_REF] is only available with light fields with narrow baseline), and the size of the datasets is sufficient for training a deep neural network for light field disparity/depth estimation. Our datasets will be available online upon paper acceptance.

Training

Training stereo pairs are randomly selected among pairs of LF views in the same row or the same column. For dense LFs, the angular distance between these two stereo views is between 2 and 8 view indices. The same distance is chosen between 1 and 3 for sparse light fields. In order to avoid imbalanced data distribution, we have managed to keep equal the occurrence of training examples for each different distance. The model is first trained with data in SLF dataset. Then, for dense light fields, the model is further finetuned using the DLF dataset. We also include the 16 additional scenes (different from those evaluated in Table 1) from the HCI 4D light field dataset to enrich our DLF training set. Due to limited GPU memory, a batch size of 4 is adopted. Data augmentation. Chromatic augmentation (variation of contrast, color and brightness) is applied with the same parameters as suggested in [START_REF] Dosovitskiy | Flownet: Learning optical flow with convolutional networks[END_REF]. However, we do not apply any geometrical transformation, e.g., translation, rotation and scaling. In our experiments, it is observed that geometrical transformation implies interpolation errors in the transformed ground truth disparity maps, which harms the learning convergence. Learning rate schedule. To finetune FlowNet 2.0, the initial learning rate is set to 0.0001 for the first 500 epochs, whereas this learning rate is applied for the first 1200 epochs to train the MSRNet. Then, for both training tasks, the learning rate is decreased by half after every 200 epochs for better convergence.

Dense light fields

We assess our algorithm against several state-of-the-art methods [START_REF] Jeon | Accurate depth map estimation from a lenslet light field camera[END_REF][START_REF] Zhang | Robust depth estimation for light field via spinning parallelogram operator[END_REF][START_REF] Huang | Empirical bayesian light-field stereo matching by robust pseudo random field modeling[END_REF][START_REF] Jiang | Depth estimation with occlusion handling from a sparse set of light field views[END_REF] with 8 test dense LFs in the HCI dataset. The methods [START_REF] Jeon | Accurate depth map estimation from a lenslet light field camera[END_REF] and [START_REF] Zhang | Robust depth estimation for light field via spinning parallelogram operator[END_REF] make use of the full 4D LF. The method [START_REF] Huang | Empirical bayesian light-field stereo matching by robust pseudo random field modeling[END_REF] considers a sparse subset of 3×3 views including the center view, whereas [START_REF] Jiang | Depth estimation with occlusion handling from a sparse set of light field views[END_REF] and our method only exploit the 4 corner views.

Jeon et al. [START_REF] Jeon | Accurate depth map estimation from a lenslet light field camera[END_REF] Zhang et al. [START_REF] Zhang | Robust depth estimation for light field via spinning parallelogram operator[END_REF] Jiang et al. [START_REF] Jiang | Depth estimation with occlusion handling from a sparse set of light field views[END_REF] Huang et al. [START_REF] Huang | Empirical bayesian light-field stereo matching by robust pseudo random field modeling[END_REF] Ours GT Image We use the same evaluation metrics defined in [START_REF] Honauer | A dataset and evaluation methodology for depth estimation on 4d light fields[END_REF][START_REF] Johannsen | A taxonomy and evaluation of dense light field depth estimation algorithms[END_REF]. MSE is mean-square-error, which penalizes large disparity errors on the object boundary, whereas BadPix(α) (the percentage of pixels having an error superior to α, α being set to small values) and Q25 (the error value *100 at the 25th percentile of the disparity estimates) measure the sub-pixel accuracy. Table 1 shows that in terms of MSE, our method is on par with [START_REF] Zhang | Robust depth estimation for light field via spinning parallelogram operator[END_REF] and [START_REF] Huang | Empirical bayesian light-field stereo matching by robust pseudo random field modeling[END_REF], and better than other reference methods. In terms of BadPix(0.01), BadPix(0.03) and Q25, our method outperforms all the reference methods by a large margin. Moreover, unlike [START_REF] Jeon | Accurate depth map estimation from a lenslet light field camera[END_REF][START_REF] Zhang | Robust depth estimation for light field via spinning parallelogram operator[END_REF][START_REF] Huang | Empirical bayesian light-field stereo matching by robust pseudo random field modeling[END_REF], the fact that our algorithm generates one disparity map per viewpoint, without exploiting the color information of the views other than the four corner views, is especially interesting for applications such as light field view synthesis.

Sparse light fields

Our algorithm has been also evaluated on 4 sparse LFs in our SLF dataset. An angular resolution of 3 × 3 is considered. In this case, methods [START_REF] Jeon | Accurate depth map estimation from a lenslet light field camera[END_REF] and [START_REF] Zhang | Robust depth estimation for light field via spinning parallelogram operator[END_REF] are no longer relevant for comparison, since they only rely on densely sampled views, and their performance drops drastically when the baseline increases. In Table 2, our algorithm performs significantly better than [START_REF] Jiang | Depth estimation with occlusion handling from a sparse set of light field views[END_REF] and [START_REF] Huang | Empirical bayesian light-field stereo matching by robust pseudo random field modeling[END_REF] in subpixel accuracy (BadPix and Q25), whereas the method [START_REF] Huang | Empirical bayesian light-field stereo matching by robust pseudo random field modeling[END_REF] excels in MSE. Note that the loss function (Eq. ( 8)) computes SAD instead of MSE. For the test LFs, our algorithm obtains better SAD measures against [START_REF] Huang | Empirical bayesian light-field stereo matching by robust pseudo random field modeling[END_REF].

CONCLUSIONS

In this paper, we have proposed a learning-based depth estimation solution both for densely and sparsely sampled light field data. The experiments show that our algorithm outperforms state-of-the-art algorithms by a large margin in most of the metrics. Our algorithm can be also naturally integrated into a light field view synthesis pipeline, since it is able to infer disparity information for a view that the color information is unknown.

Fig. 1 :

 1 Fig. 1: Disparity estimation for LF anchor views by using the finetuned FlowNet 2.0. Vertical pairs are rotated 90 degrees to become horizontal pairs. An example is given on the right for the view I1,1.

Fig. 2 :

 2 Fig. 2: Network structure for light field multi-view disparity refinement. Given the coarse disparity maps d r and d r computed using image pairs, the fusion of these maps is performed by exploiting the warping errors Eocc and Enocc. The refinement of the fused disparity map Dr is learned with multi-scale residual learning.

  d r (p) and d r (p) by minimizing the error term Enocc or Eocc, respectively: D nocc r (p) = argmin d r (p),d r (p) Enocc(d r , p), Enocc(d r , p) , argmin d r (p),d r (p) Eocc(d r , p), Eocc(d r , p) .

  B = |D • C(Dr) -D • C(Dr) D • C(Ir)|. D and C are respectively dilation and canny edge detection operators. The symbol • denotes function composition and is for Hadamard product.

  l 2 -norm L 2 (D, D gt ) = p G(p) 2 withG(p) = ∇xD(p) -∇xD gt (p), ∇yD(p) -∇yD gt (p) .[START_REF] Ting-Chun | Occlusion-aware depth estimation using light-field cameras[END_REF]

Fig. 3 :

 3 Fig. 3: Visual comparison of the estimated disparity maps on center view.

Table 1 :

 1 Quality evaluation of the estimated disparity maps on center view for dense LFs. The best results are marked in bold. MonasRoom 0.76 0.58 0.53 0.56 0.38 46.0 42.5 38.6 45.5 25.2 22.1 17.8 18.6 17.8 11.4 0.38 0.34 0.33 0.35 0.24 Sideboard 4.16 0.92 10.31 1.31 1.60 67.8 64.3 61.7 79.8 48.8 39.3 31.0 37.5 46.4 23.2 0.74 0.66 1.26 0.51 0.37

		MSE*100	BadPix(0.01)	BadPix(0.03)	Q25
	Light fields	[2]

[START_REF] Zhang | Robust depth estimation for light field via spinning parallelogram operator[END_REF]

[START_REF] Jiang | Depth estimation with occlusion handling from a sparse set of light field views[END_REF]

[10] Ours

[START_REF] Jeon | Accurate depth map estimation from a lenslet light field camera[END_REF]

[5]

[START_REF] Jiang | Depth estimation with occlusion handling from a sparse set of light field views[END_REF]

[10] Ours

[START_REF] Jeon | Accurate depth map estimation from a lenslet light field camera[END_REF]

[5]

[START_REF] Jiang | Depth estimation with occlusion handling from a sparse set of light field views[END_REF]

[10] Ours

[START_REF] Jeon | Accurate depth map estimation from a lenslet light field camera[END_REF]

[5] [11] [10] Ours StillLife 2.02 1.72 2.56 1.16 1.14 81.2 76.2 71.3 74.4 71.5 51.0 32.1 25.0 37.1 24.5 1.36 1.02 0.87 0.86 0.88 Buddha 1.13 0.97 0.82 0.40 0.46 57.7 41.2 34.9 51.3 25.8 24.4 14.8 12.3 13.4 6.6 0.51 0.34 0.31 0.52 0.28 Butterfly 4.79 0.74 1.84 0.70 0.54 82.5 78.9 70.8 82.4 62.9 49.1 48.5 36.0 50.8 28.7 1.47 1.22 0.85 1.28 0.66 Boxes 14.15 8.23 12.71 10.05 12.48 72.7 62.3 65.8 83.6 60.5 45.5 28.1 37.7 57.1 32.8 0.89 0.62 0.68 1.54 0.55 Cotton 9.98 1.44 1.18 1.23 0.67 60.5 41.7 42.6 72.1 29.6 23.3 11.1 10.7 33.7 8.0 0.59 0.36 0.42 0.89 0.25 Dino 1.23 0.29 0.88 0.53 0.50 76.6 57.5 49.1 80.9 35.9 48.4 17.9 20.0 48.0 12.6 1.08 0.55 1.32 0.42 0.29 Average 4.78 1.86 3.85 1.99 2.22 68.1 58.1 54.4 71.2 45.0 37.9 25.2 24.7 38.0 18.5 0.88 0.64 0.62 0.80 0.44

Table 2 :

 2 Quality evaluation of the estimated disparity maps on center view for sparse LFs.

		MSE	BadPix(0.1)	Q25
	Light fields [11] [10] Ours [11] [10] Ours [11] [10] Ours
	Furniture	1.94 0.38 0.78 41.3 61.3 22.0 2.52 6.17 1.10
	Lion	0.87 0.08 0.15 59.5 21.4 8.0 4.47 2.51 0.61
	Toy bricks 1.10 0.18 0.44 44.6 36.0 16.6 3.61 2.72 0.94
	Electro devices 0.63 0.18 0.23 43.4 55.5 24.5 2.71 4.93 1.35
	Average	1.14 0.21 0.40 47.2 43.6 17.8 3.33 4.08 1.00

& ↓ 2 ( ↓