Fatma Hawary

Guillaume Boisson

Christine Guillemot

Philippe Guillotel

COMPRESSIVE 4D LIGHT FIELD RECONSTRUCTION USING ORTHOGONAL FREQUENCY SELECTION

Keywords: Light field, reconstruction, non-regular sampling, sparsity

We present a new method for reconstructing a 4D light field from a random set of measurements. A 4D light field block can be represented by a sparse model in the Fourier domain. As such, the proposed algorithm reconstructs the light field, block by block, by selecting frequencies of the model that best fits the available samples, while enforcing orthogonality with the approximation residue. The method achieves a very high reconstruction quality, in terms of Peak Signal-to-Noise Ratio (PSNR). Experiments on several datasets show significant quality improvements of more than 1dB compared to state-of-the-art algorithms.

INTRODUCTION

Light field imaging has recently gained in popularity due to the emergence of acquisition devices and the potential for a variety of computational photography and computer vision applications. Nevertheless, existing acquisition devices either provide a low-resolution light field by multiplexing several views on a single 2D sensor image, or require capturing multiple images to generate a high-resolution light field. This limitation goes against industrial needs for increasing image resolutions. In parallel, in the field of compressive sensing [START_REF] Candès | Sparse Signal and Image Recovery from Compressive Samples[END_REF], it has been shown that it is possible to recover a signal from only a few samples (or measurements) provided this signal is sparse in a particular domain. Light field camera architectures have hence been proposed based on compressive acquisition schemes with the goal of overcoming the spatio-angular trade-off of plenoptic cameras.

Programmable aperture approaches exploit the fast multiple-exposure feature of digital sensors to sequentially capture subsets of light rays [START_REF] Liang | Programmable aperture photography: multiplexed light field acquisition[END_REF]. The authors in [START_REF] Ashok | Compressive Light Field Imaging[END_REF] propose instead two compressive acquisition architectures, one exploiting spatial correlations and another exploiting angular correlations. In [START_REF] Xu | A high-resolution light field camera with dual-mask design[END_REF], two attenuation masks are used, one placed at the aperture and the other one in front of the 2D photosensor. Similarly, a random-coded mask is used in [START_REF] Babacan | Compressive Light Field Sensing[END_REF] to capture random linear combinations of angular samples. An architec-ture using coded projections on the sensor image is described in [START_REF] Marwah | Compressive Light Field Photography using Overcomplete Dictionaries and Optimized Projections[END_REF], where the light field is reconstructed using sparse methods with an over-complete dictionary learned with K-SVD. Casting the light field acquisition problem in a compressive sensing framework requires learning dictionaries well suited for the data at hand. The dictionary learning problem for light field compressive sensing has been addressed in [START_REF] Miandji | Compressive Image Reconstruction in Reduced Union of Subspaces[END_REF]: the authors train an ensemble of separable 2D dictionaries corresponding to a reduced union of subspaces (RUS) in which the input data is supposed to reside. Alternatively, a perspective shifting based dictionary is defined in [START_REF] Chen | Light Field Compressed Sensing Over a Disparity-Aware Dictionary[END_REF] to sparsely represent light fields and two disparity estimation methods are proposed for the reconstruction of the light field content.

In this paper, we address the problem of recovering a 4D light field signal from a subset of random measurements. The random sampling avoids making any further hypotheses. We exploit the assumption that the light field data is sparse in the Fourier domain [START_REF] Hawary | Scalable Light Field Compression Scheme Using Sparse Reconstruction and Restoration[END_REF], meaning that the Fourier transform of the light field can be expressed as a linear combination of a small number of Fourier basis functions. The reconstruction algorithm therefore searches for these bases (i.e. their frequencies) which best represent the 4D Fourier spectrum of the sampled light field. The method is inspired from the Frequency Selective Reconstruction (FSR) approach in [START_REF] Seiler | Resampling Images to a Regular Grid From a Non-Regular Subset of Pixel Positions Using Frequency Selective Reconstruction[END_REF], but extended to 4D light fields. The method is further improved by introducing an orthogonality constraint on the residue in the same spirit of Orthogonal Matching Pursuit (OMP) [START_REF] Cai | Orthogonal matching pursuit for sparse signal recovery[END_REF] compared to Matching Pursuit. The experimental results on both synthetic and real contents show that the proposed method provides a high-quality reconstruction, even with a very low number of input samples. The improved version of the algorithm achieves gains of up to 1.9dB in PSNR compared to 4D-FSR, and outperforms the state-of-the-art method in [START_REF] Miandji | Compressive Image Reconstruction in Reduced Union of Subspaces[END_REF].

LIGHT FIELD RESAMPLING

Problem statement

We consider a 4D randomly-sampled light field, and we aim at estimating the 4D signal that best fits the available samples. The light field reconstruction is conducted per 4D hyper- The imaginary unit:

j 2 = -1 z * = (z) -j. (z) Complex conjugate of z i Iteration index Ω = 1; K × 1; L × 1; M × 1; N Local light field domain P = |Ω| = K.L.M.N Number of samples in Ω A ⊂ Ω Subset of known samples B ⊂ Ω Subset of unknown samples C ⊂ Ω Subset of reconstructed samples p = (k, l, m, n) ∈ Ω A pixel's position within Ω f : Ω → R Local light field g (i) : Ω → R Approximation model at iteration i w : Ω → R Weighting function r (i) : Ω → R Weighted residue at iteration i ϑ = (s, t, u, v)
A frequency in the 4D spectrum

ϕ ϑ : Ω → R A 4D Fourier basis function Θ (i) = {ϑ 1 , . . . , ϑ i } Frequency subset at iteration i X ϑ = p x[p]ϕ * ϑ [p]
Capitals denote Fourier transforms block, in the order of decreasing density of known samples. Each hyper-block to reconstruct is regarded as the core of a surrounding area in the spatial and angular dimensions (see Fig. 1). Let Ω be the 4D domain made of the hyper-block and its local surroundings: Ω spans over M × N views with colocated K × L patches. Let A, B and C respectively denote the local subsets of the known, unknown, and reconstructed samples: 1 summarizes the notations used throughout the paper.

Ω = A ∪ B ∪ C. Table
A weighting function w is introduced to discriminate known and reconstructed samples from unknown samples in this area:

w[p] = w[k, l, m, n] =        ρ √ k 2 +l 2 s ρ √ m 2 +n 2 a for p ∈ A 0 for p ∈ B δ.ρ √ k 2 +l 2 s ρ √ m 2 +n 2 a for p ∈ C (1)
where

k = k -K-1 2 , l = l -L-1 2 , m = m -M -1 2 , n = n -N -1
2 , and 0 < ρ s , ρ a , δ < 1. The factors ρ s and ρ a respectively determine the weighting decay in the spatial dimension (within a view) and in the angular (cross-view) dimension. The parameter δ is used to differentiate the contribution of reconstructed data from the available data.

Eventually, let f be the local light field defined on Ω. The algorithm generates its sparse approximation model

g[p] = ϑ∈Θ c ϑ ϕ ϑ [p] (2)
as a weighted combination of Fourier basis functions defined as:

ϕ ϑ [p] = ϕ stuv [k, l, m, n] = e 2πj( ks K + lt L + mu M + nv N ) (3)

Iterative light field reconstruction by frequency selection

Let r (i) be the weighted residue of the approximation model with respect to the signal at iteration i: A global overview of the reconstruction method is provided in Algorithm 1. The model is generated iteratively by adding a new basis function at each iteration. Two main steps are conducted at each iteration: the selection of the basis function that minimizes the residual error, followed by the update of the model and the residue.

r (i) = (f -g (i) ).w (4)
The basic idea of the algorithm FSR [START_REF] Seiler | Resampling Images to a Regular Grid From a Non-Regular Subset of Pixel Positions Using Frequency Selective Reconstruction[END_REF] is similar to Matching Pursuit [START_REF] Mallat | Matching Pursuits with Time-frequency Dictionaries[END_REF] in the sense that it aims at sparsely approximating a signal by finding the best matching projections of the data onto the span of a Fourier basis. As such, the selected basis function at iteration i is the one on which the projection of the residue r (i-1) is maximal (line 5. in Algorithm 1). The corresponding expansion coefficient c

(i)
ϑi is computed by minimizing the weighted residual energy E (i)

w :

E (i) w = p∈Ω (f [p] -g (i-1) [p] -c (i) ϑi ϕ ϑi [p]) 2 .w[p] (5)
According to the authors of [START_REF] Seiler | Resampling Images to a Regular Grid From a Non-Regular Subset of Pixel Positions Using Frequency Selective Reconstruction[END_REF], minimizing

E (i)
w with respect to c (i) ϑi and c (i) * ϑi leads to the following equation:

c (i) ϑ i = p∈Ω r (i-1) [p]ϕ * ϑ i [p] p∈Ω w[p] (6)
Once the basis function is selected, the sparse model as well as the residue are updated (lines 7. and 8. in Algorithm 1). The algorithm then proceeds to the next iteration where a new basis function is selected, and so on, until a predefined number of iterations is reached. Since the basis functions are not orthogonal if evaluated with the weighting function, a constant factor γ is introduced in [START_REF] Seiler | Resampling Images to a Regular Grid From a Non-Regular Subset of Pixel Positions Using Frequency Selective Reconstruction[END_REF] in each expansion coefficient, to compensate the orthogonality deficiency. So far, the minimization criterion takes into account the residue projection only on the current selected basis function, and thus, does not ensure having the best approximation with all basis functions selected so far. Consequently, a basis function can Initialization: g (0) = 0; r (0) = f.w 4:

for iteration i ← 1, max iterations do 5:

ϑ i = argmax ϑ ϕ ϑ , r (i-1)

6:

Compute c (i) ϑi 7:

g (i) = g (i-1) + c (i)
ϑi ϕ ϑi 8:

r (i) = r (i-1) -c (i) ϑi ϕ ϑi w 9:
Update the light field on the hyper-block be selected more than once, and a high number of iterations (100 in [START_REF] Seiler | Resampling Images to a Regular Grid From a Non-Regular Subset of Pixel Positions Using Frequency Selective Reconstruction[END_REF]) is required to reach a fair reconstruction quality.

To overcome these limitations, we introduce a better criterion for approximating the residual error. By incorporating all the selected basis functions in the update of the approximated residue, a more accurate reconstruction can be obtained with a reduced number of iterations. This idea has been discussed in the context of image coding in [START_REF] Kaup | Coding of segmented images using shape-independent basis functions[END_REF], but has never been analyzed, nor used, to sparsely reconstruct randomly sampled light fields.

Hermitian symmetry and orthogonality criterion

Since the light field signal is Hermitian, we explicitly gather the conjugate contributions of opposite frequencies at each iteration. We furthermore take into account the impact of the new selected frequency upon every frequency selected so far:

g (i) = g (i-1) + ϑ∈Θ (i) 1 2 ∆c ϑ ϕ ϑ + ∆c * ϑ ϕ * ϑ (7)
The basis function selection is conducted in the same way as in the previous section. However, the new minimization criterion with respect to each ∆c ϑ yields a system of equations whose solution ensures that the residue is orthogonal to each basis function selected so far:

p∈Ω w[p].ϕ ϑ [p] ϑ ∈Θ (i) 1 2 ∆c ϑ ϕ ϑ [p] + ∆c * ϑ ϕ * ϑ [p] = p∈Ω r (i-1) [p]ϕ ϑ [p], ∀ϑ ∈ Θ (i) (8)
Equivalently, minimizing with respect to ∆c * ϑ leads to the complex conjugate of the system in Eq.8. Thanks to this new criterion, a basis function is selected only once, and the model g is completely recalculated at each iteration.

In the Fourier domain, the previous equation can be expressed as:

ϑ ∈Θ (i) 1 2 ∆c ϑ W * ϑ +ϑ + ∆c * ϑ W ϑ -ϑ = R (i-1) * ϑ , ∀ϑ ∈ Θ (i) (9)
where W and R refer to the Fourier transforms of w and r, respectively. The objective is to estimate updates of all the expansion coefficients {∆c (i) ϑt } t=1..i , i.e. to resolve the following equation:

∆c (i) = 2 W (i) -1 R (i-1) (10)
where we define the (2i -1) vectors ∆c (i) and R (i-1) , and the matrix W (i) of size (2i -1) × (2i -1) as follows:

∆c (i) =             (∆c (i) ϑ 1 ) . . . (∆c (i) ϑ i ) (∆c (i) ϑ 2 ) . . . (∆c (i) ϑ i )             R (i-1) =             (R (i-1) ϑ 1 ) . . . (R (i-1) ϑ i ) (R (i-1) ϑ 2 ) . . . (R (i-1) ϑ i )             (11)
W (i) = W (i) 11
W (i) 12
W (i) 21 W (i) 22
with:

W (i) 11 = (W ϑx+ϑy + W ϑx-ϑy ) (x,y)∈ 1;i × 1;i W (i) 12 = (W ϑx+ϑy -W ϑx-ϑy ) (x,y)∈ 1;i × 2;i W (i) 21 = (W ϑx+ϑy + W ϑx-ϑy ) (x,y)∈ 2;i × 1;i W (i) 22 = (W ϑx-ϑy -W ϑx+ϑy ) (x,y)∈ 2;i × 2;i
Eventually, the parametric model and the residue are updated in the Fourier domain according to all the {1..i} selected basis functions:

G (i) ϑ = G (i-1) ϑ + 1 2 P.∆c (i) ϑ G (i) -ϑ = G (i-1) -ϑ + 1 2 P.∆c * (i) ϑ , ∀ϑ ∈ Θ (i) (12)
R (i) ϑ = R (i-1) ϑ - ϑ ∈Θ (i) 1 2 ∆c ϑ W * ϑ -ϑ + ∆c * ϑ W ϑ +ϑ , ∀ϑ (13)

EXPERIMENTAL RESULTS

We show experimental results of the reconstruction of subsampled light fields by the iterative 4D-FSR algorithm presented beforehand, using the same parameter values as in [START_REF] Seiler | Resampling Images to a Regular Grid From a Non-Regular Subset of Pixel Positions Using Frequency Selective Reconstruction[END_REF], with iteration number fixed at 100, as well as by the Orthogonal Frequency Selection (OFS) introduced in section 2.3, but with a number of iterations reduced to 30. The parameters ρ s , ρ a and δ are respectively set to 0.7, 0.5 and 0.5. The reconstruction quality is evaluated in terms of PSNR averaged on all views for different sampling rates. StillLife and Butterfly [START_REF] Wanner | Datasets and benchmarks for densely sampled 4d light fields[END_REF]. The proposed method 4D-OFS achieves a high-quality reconstruction of the tested contents, with PSNR reaching 40dB for some contents, and significant average gains ranging from 0.52dB to 1.9dB in comparison to 4D-FSR. Fig. 3 compares the reconstruction quality of a view from the light field StillLife [START_REF] Wanner | Datasets and benchmarks for densely sampled 4d light fields[END_REF], reconstructed by both methods at a sampling rate of 4%. One can see that the 4D-FSR presents more artifacts on the edges. Additionally, Fig. 4 compares these methods with the state-of-the-art compressive sensing algorithm (RUS) [START_REF] Miandji | Compressive Image Reconstruction in Reduced Union of Subspaces[END_REF] using non-overlapping blocks, on the light field Dragon [START_REF] Marwah | Compressive Light Field Photography using Overcomplete Dictionaries and Optimized Projections[END_REF]. 4D-OFS achieves a significant gain in reconstruction quality compared to RUS, especially for a low number of samples (2.7dB at 4%). A comparison of the visual quality of the reconstruction can be seen in Fig. 5.

CONCLUSION

We introduced a new iterative block-wise algorithm to reconstruct randomly sampled light fields. The work falls within the effort of addressing the challenge of dense light field acquisition with compressive sensing schemes. The method selects the Fourier basis functions that best approximate the available samples, and generates a sparse model of the signal. We then improve the method by presenting an orthogonal version which ensures that the approximation residue lies in a subspace orthogonal to the one spanned by the already selected basis functions. Experimental results show that the method achieves a high-quality reconstruction, and demonstrate the advantage of the improved version, while using a reduced number of iterations. A comparison with a learningbased state-of-the-art method demonstrates PSNR gains, even for low sampling rates. Future work will focus on extrapolating the approximation into the continuous Fourier domain.

Fig. 1 :

 1 Fig. 1: 4D hyper-block (outlined in red) and its local spatio-angular neighborhood (K = L = 6 and M = N = 5).

Algorithm 1 field reconstruction 1 :

 11 Light Split the light field into 4D hyper-blocks 2: for each hyper-block in the decreasing density order do 3:

Fig. 2 :Fig. 3 :

 23 Fig. 2: Tested light fields. From left to right: Flower1, Cars, Crystal, Lego Truck, StillLife, Butterfly.

Fig. 4 :Fig. 5 :

 45 Fig. 4: PSNR comparison of 4D-FSR, 4D-OFS methods and RUS [7] for different sampling rates.

Table 1 :

 1 Table of notations

	j

  Table 2 represents results of tests carried out on various light fields: plenoptic contents Flower1 and Cars from [14], multi-view contents Crystal and Lego Truck 1 as well as synthetic datasets

Table 2 :

 2 PSNR in dB from light fields reconstruction by the proposed method 4D-OFS compared to FSR[START_REF] Seiler | Resampling Images to a Regular Grid From a Non-Regular Subset of Pixel Positions Using Frequency Selective Reconstruction[END_REF] extended to 4D. .44 30.47 31.67 34.44 36.38 36.15 37.17 28.20 29.02 39.59 40.50 10% 32.77 33.31 30.91 32.11 34.84 36.77 36.57 37.62 28.74 29.77 39.83 40.87 20% 33.34 33.82 31.47 32.64 35.48 37.38 37.21 38.26 29.66 30.69 40.46 41.41 40% 34.78 35.31 32.73 33.88 36.71 38.57 38.46 39.49 31.33 32.39 41.81 43.29

			Plenoptic data			Multi-view data			Synthetic data
	Sampling rate	Flower1	Cars	Crystal	Lego Truck	StillLife	Butterfly
		FSR	OFS	FSR	OFS	FSR	OFS	FSR	OFS	FSR	OFS	FSR	OFS
	4%	31.92 32

http://lightfield.stanford.edu/lfs.html