Local Fourier Slice Photography

TL;DR

This paper introduces local Fourier slice photography, enabling refocused image reconstruction directly from sparse wavelet representations of light fields, significantly reducing storage and computation costs while maintaining image quality.

Contribution

It presents a novel wavelet-based method that allows exact refocused image reconstruction from sparse light field data without dense pixel reconstruction.

Findings

Reduces storage requirements for light field images

Enables fast, exact refocused image reconstruction from sparse data

Scales linearly with visual complexity and region size

Abstract

Light field cameras provide intriguing possibilities, such as post-capture refocus or the ability to synthesize images from novel viewpoints. This comes, however, at the price of significant storage requirements. Compression techniques can be used to reduce these but refocusing and reconstruction require so far again a dense pixel representation. To avoid this, we introduce local Fourier slice photography that allows for refocused image reconstruction directly from a sparse wavelet representation of a light field, either to obtain an image or a compressed representation of it. The result is made possible by wavelets that respect the "slicing's" intrinsic structure and enable us to derive exact reconstruction filters for the refocused image in closed form. Image reconstruction then amounts to applying these filters to the light field's wavelet coefficients, and hence no reconstruction of a dense pixel representation is required. We demonstrate that this substantially reduces storage requirements and also computation times. We furthermore analyze the computational complexity of our algorithm and show that it scales linearly with the size of the reconstructed region and the non-negligible wavelet coefficients, i.e. with the visual complexity.

Equations

$$\displaystyle\hat{\psi}_{s}({\xi})\equiv\hat{\psi}_{jkt}({\xi})=\frac{2^{j}}{2\pi}\Big{(}\sum_{n}\beta_{j,n}^{t}\,e^{in\theta_{{\xi}}}\Big{)}\,\hat{h}(2^{-j}|{\xi}|)\,e^{-i\langle{\xi},2^{j}{k}\rangle}$$

$${\psi }_s(x)\equiv {\psi }_{jkt}(x)=\frac{2^j}{2\pi }\sum _n{i}^n{\beta }_{j,n}^t{e}^{in{\theta }_x}{h}_n(|2^{j}x-k|)$$

$$f(x)$$

$$\displaystyle I(x,y)=\frac{1}{\alpha^{2}}\!\int_{\mathbb{R}_{u}}\!\int_{\mathbb{R}_{v}}\!\!\ell\big{(}x/\alpha+(1\!-\!1/\alpha)u,y/\alpha+(1\!-\!1/\alpha)v,u,v\big{)}\,du\,dv$$

$$\displaystyle I(x,y)=\int_{\mathbb{R}_{u}}\!\int_{\mathbb{R}_{v}}\ell\big{(}S_{\alpha}(x,u)^{T};S_{\alpha}(y,v)^{T}\big{)}\,du\,dv$$

$$\displaystyle S_{\alpha}=\left(\begin{array}[]{cc}1/\alpha&1-1/\alpha\\ 0&1\end{array}\right).$$

$$\ell (x,u;y,v)=\sum _{(s,r)\in \mathcal{L}}{\ell }_{sr}{\psi }_s(x,u){\psi }_r(y,v)$$

$$I(x,y)$$

$$\displaystyle=\sum_{(s,r)\in\mathcal{L}}\ell_{sr}\int_{\mathbb{R}_{u}}\!\psi_{s}\big{(}S_{\alpha}(x,u)^{T}\big{)}\,du\int_{\mathbb{R}_{v}}\!\psi_{r}\big{(}S_{\alpha}(y,v)^{T}\big{)}\,dv.$$

$$\displaystyle g(x)=\int_{\mathbb{R}_{u}}f\big{(}S_{\alpha}(x,u)^{T}\big{)}\,du$$

$$\displaystyle g(x)=\sum_{s\in\mathcal{S}}f_{s}\,\int_{\mathbb{R}_{u}}\psi_{s}\big{(}S_{\alpha}(x,u)^{T}\big{)}\,du.$$

$$\displaystyle\int_{\mathbb{R}_{u}}\!\psi_{s}\big{(}S_{\alpha}(x,u)^{T}\big{)}\,du=\mathcal{F}_{x}^{-1}\left(\alpha^{-1}\,\hat{\psi}_{s}\big{(}S_{\alpha}^{-T}(\xi_{x},\xi_{u})^{T}\big{|}_{\xi_{u}=0}\big{)}\right)$$

$${\int }_{{\mathbb{R}}_u}$$

$$\displaystyle\times\underbrace{\int_{\mathbb{R}_{\xi_{x}}}e^{in\theta_{S_{\alpha}^{-T}\xi_{0}}}\,\hat{h}\big{(}2^{-j}|S_{\alpha}^{-T}{\xi}_{0}|\big{)}\,e^{-i\langle{\xi_{0}},2^{j}S^{-1}{k}\rangle}\,e^{i\xi_{x}x}\,d\xi_{x}}_{\displaystyle\zeta_{s}^{\alpha,n}(x)=\zeta_{j_{s},t_{s}}^{\alpha,n}\big{(}x-P_{x}(2^{j}S^{-1}{k_{s}})\big{)}},$$

$${\zeta }_s^{\alpha }(x)\equiv \frac{2^j{\alpha }^{-1}}{(2\pi {)}^{3/2}}\sum _n{\beta }_{j,n}^t{\zeta }_s^{\alpha ,n}(x)={\int }_{{\mathbb{R}}_u}$$

$$\displaystyle g(x)=\sum_{s\in\mathcal{S}}f_{s}\,\int_{\mathbb{R}_{u}}\psi_{s}\big{(}S_{\alpha}(x,u)^{T}\big{)}\,du=\sum_{s\in\mathcal{S}}f_{s}\,\zeta_{s}^{\alpha}(x).$$

$$g(x)$$

$$I(x,y)$$

$${\ell }_{qp}$$

$$=\sum _{(s,r)\in \mathcal{L}}{\ell }_{sr}{⟨{\zeta }_s^{\alpha }(x),{\psi }_q^1(x)⟩}_x{⟨{\zeta }_r^{\alpha }(y),{\psi }_p^1(y)⟩}_y.$$

$${\gamma }_{sq}=⟨{\zeta }_s^{\alpha }(x),{\psi }_q^1(x)⟩$$

$${\ell }_{qp}$$

$$K_{sr}=2c_{\zeta }{2}^{-j_s-j_r}{W}^{2}R$$

$$\displaystyle\big{|}\mathcal{L}_{j_{s},j_{r}}\big{|}=4\cdot\,2^{-2(j_{max}-j_{s})}\,2^{-2(j_{max}-j_{r})}\,N_{x}^{2}\,N_{u}^{2}\,T_{j_{s}}\,T_{j_{r}}\,|A|$$

$$K_{j_s,j_r}$$

$$K_{j_{max}}$$

$$K_{j_{max}}^{\alpha }$$

$${\mathrm{cr}}_{j_s,j_r}^{ϵ}=\frac{|{\mathcal{L}}_{j_s,j_r}|}{|{\mathcal{L}}_{j_s,j_r}^{\alpha ,ϵ}|}.$$

$$K_{j_{max}}^{\alpha ,ϵ}$$

$$\displaystyle K_{{j_{max}}}^{\alpha,\epsilon}=\mathcal{O}\Big{(}\,P_{j_{max}}\,\frac{N_{x}^{2}\,N_{u}^{2}}{cr^{\epsilon}}\,|A|\,R\Big{)}.$$

$${\zeta }_s^{\alpha ,n}(x)$$

$$\times e^{-in\left(\pi -{\tan }^{-1}\left(\left(1-\frac{1}{\alpha }\right)\alpha \right)\right)}{\left({\alpha }^2{x}^2\right)}^{-\frac{i\pi }{\log (4)}}$$

Full text

light field camera, Fourier slice theorem, wavelets

††journalyear: 2019††copyright: acmlicensed††conference: Siggraph ’19; ; ††booktitle: Siggraph ’19††price: 15.00††doi: 10.1145/1122445.1122456††isbn: 978-1-4503-9999-9/18/06††ccs: Computing methodologies Computational photography††ccs: Computing methodologies Image compression††ccs: Mathematics of computing Computation of transforms

1. Introduction

Light field cameras, which record the full four-dimensional plenoptic function, open up many possibilities for both consumer, e.g. (Ng et al., 2005), and professional applications, e.g. (Levoy et al., 2006). Prime examples are post-capture refocus and the ability to obtain images where every depth is in focus. The possibilities, however, come at the price of considerable storage requirements for the light field data sets. Compression techniques can alleviate these but image reconstruction and light field processing typically require again a dense representation.

To avoid this, we propose local Fourier slice photography, an algorithm to compute refocused images directly from a light field’s sparse wavelet representation. Our work draws inspiration from Ng’s seminal Fourier slice photography (2005) where image reconstruction is performed in the frequency domain using the projection slice theorem. We combine this work with a recent advancement of the slice theorem (Lessig, 2018a) that uses carefully chosen wavelets to allow for an efficient projection from a signal’s compressed wavelet representation. To apply this result to refocused image reconstruction, we extend it to a sheared, local projection slice equation that establishes closed-form, shear-dependent reconstruction kernels for the projected signal. With these, a refocused image can be obtained directly from a light field’s compressed wavelet coefficients using an inverse transform. We also derive an extension that enables one to directly obtain sparse, refocused images from a sparse light field data set. Our experimental results confirm that our approach yields high fidelity, refocused images directly from compressed light fields without the need to obtain a dense pixel representation. They also demonstrate that errors that arise at high compression rates mainly manifest themselves through lost high-frequency detail, i.e. without distracting artifacts.

The sparsity that reduces storage requirements also reduces computational costs. We show this experimentally and verify it through a theoretical analysis that establishes a linear dependence on the number of nonzero wavelet coefficients. Because of the spatial localization of the wavelets, the costs of our technique depend on the light field’s angular resolution. This was not the case for Fourier slice photography (Ng, 2005) although there one cannot easily take advantage of redundancy in the data. The localization also enables us to obtain all-focus images, which is not possible using Fourier slice photography.

A conceptual overview of our approach is provided in Fig. 2 and the computations required in an implementation are summarized in Algorithm 1. The remainder of the paper, which expounds on the details, is structured as follows. After reviewing related work in the next section, we provide in Sec. 3 the necessary background on the polar wavelets that are used in our work. In Sec. 4 we derive the sheared local Fourier slice equation and develop our technique to obtain a refocused image directly from a light field’s sparse wavelet representation. Experimental results on refocused image reconstruction and all-focus images as well as details on our reference implementation are presented in Sec. 5. We conclude in Sec. 7 with a discussion of possible directions for future work.

2. Related Work

In computer graphics, light fields were introduced by Levoy and Hanrahan (1996) and Gortler et al. (Gortler et al., 1996). Initially, these were mainly of academic interest, e.g. (Chai et al., 2000), but in the 2000s practical means to capture real-world light field data sets became available (Wilburn et al., 2005; Ng et al., 2005; Venkataraman et al., 2013). With these, the generation, processing, and display of light fields has become an important research direction (Ihrke et al., 2016; Wu et al., 2017). In the following, we will therefore focus on related work most pertinent to our own.

Ng (2005) showed that post-capture refocus can be formulated in the frequency domain using the Fourier projection slice theorem and that this provides an asymptotic speedup compared to the pixel domain. Our work is inspired by Ng’s and we combine it with a recent result in optics (Lessig, 2018a) that extends the slice theorem to a spatially localized form using wavelets. The extension relies on the use of polar wavelets, which is a family of wavelets defined separably in polar coordinates in the Fourier domain (Unser and Van De Ville, 2010; Unser et al., 2011; Unser and Chenouard, 2013). These include a wide range of steerable wavelets (Perona, 1991; Freeman and Adelson, 1991; Simoncelli and Freeman, 1995) as well as curvelets (Candès and Donoho, 2005a, b) and ridgelets (Candès and Donoho, 1999; Donoho, 2000). The present work also exploits the separability of polar wavelets in polar frequency coordinates and we extend (Lessig, 2018a) to include the shearing that implements refocusing. We also benefit from the sparsity available with curvelet-like constructions (Candès and Donoho, 2004; Donoho, 2000), which ensures an efficient sparse representation of light field data sets.

Light field imaging in the Fourier domain was also considered by Shi et al. (2014). Their objective was to circumvent the sparsity degradation that results when the discrete instead of the continuous Fourier transform is used in numerical calculations. This is no issue for our technique since our reconstruction kernels are obtained in the continuous Fourier domain. Furthermore, while Shi et al. require a nonlinear optimization to obtain sparsity we use simple thresholding and rely on the compatibility of polar wavelets with the structure of natural images in frequency space (Candès and Donoho, 2005a; Mallat, 2009, Ch. 9). The design of of light field cameras and their sensors has been analyzed comprehensively by Liang and Ramamoorthi (2015). Although these authors also perform their analysis in the Fourier domain, this work is orthogonal to ours since we assume we have a preprocessed light field data set as input.

For image synthesis, the frequency representation of the light field has been analyzed in a series of papers starting with the seminal work by Durand et al. (2005). The work showed, for example, that the shearing that implements refocusing in the Fourier domain is the general expression for the transport of the light field in the two plane parametrization (Chai et al., 2000). Closest to our work are Fourier-based approaches for depth of field rendering (Soler et al., 2009; Lehtinen et al., 2011; Belcour et al., 2013). This work, however, aims at finding optimal sampling rates for the light field in Monte Carlo renderers while we assumes a (largely noise free) light field on the camera is available. Because of the curse of dimensionality, image synthesis applications also typically do not employ an explicit basis representation of the light field, which is a key part of the present work. Sen, Darabi, and Xiao (Sen et al., 2011) used compressive sensing to reduce the sampling rate for depth of field rendering. The polar wavelets employed in our work provide a sparse representation for image and light field data that would be well suited for compressive sensing.

Vagharshakyan, Bregovic, and Gotchev (2018) proposed the use of shearlets for light field reconstruction from a limited set of perspective views. Shearlets can be seen as a stereographic projection of polar wavelets with the directional localization controlled by the parabolic scaling also used for curvelets. The authors exploit the sparsity afforded by the shearlet transform to obtain an efficient algorithm for the reconstruction. However, they do not exploit that image reconstruction is naturally formulated in polar coordinates in the frequency domain, which at the heart of our work. In fact, to the best of our knowledge, with shearlets no closed form solution for the reconstruction kernels would be available.

Learning-based techniques for image reconstruction from light fields have also received considerable attention in recent work, e.g. (Levin et al., 2008; Kalantari et al., 2016; Yoon et al., 2015; Xu et al., 2018). The objective there is typically to perform the reconstruction from a reduced set of measurements, and it is hence orthogonal to our work. In the spirit of (Vagharshakyan et al., 2018), we believe that our polar wavelet representation of a light field might provide a useful pre-processing for learning-based reconstruction since it removes redundancy while respecting the intrinsic structure.

A variety of approaches for the compression of light field data sets have been proposed in the literature (Viola et al., 2017; Wu et al., 2017), for example adapting techniques used for image compression, e.g. (Alves et al., 2018); developing custom ones for slices or the full 4D light field, e.g. (Aggoun, 2006; Aggoun and Mazri, 2008; Kundu, 2012; Conti et al., 2014); or extending video compression schemes by exploiting that a light field can be seen as a sequence of images recorded from a set of nearby vantage points, e.g. (Vieira et al., 2015; Dai et al., 2015). Our use of polar wavelet for transform coding is dictated by our objective to refocus from the compressed representation. However, we do not provide a full compressions technique, e.g. we do not consider the choice of color space, gamma correction, and quantization.

3. Polar Wavelets

Wavelets are functions that are well localized in both the spatial and frequency domain. Through this, they enable an efficient and sparse representation of signals and can reduce the computational costs of numerical computations, e.g. (DeVore, 2006; Stevenson, 2009). In multiple dimensions, wavelets are typically constructed as tensor products of one dimensional ones. Polar wavelets (Unser and Chenouard, 2013), in contrast, are defined separably in polar coordinates in Fourier space, which leads to many desirable properties.

Polar wavelets are constructed using a compactly supported radial window $\hat{h}(|{\xi}|)$, which controls the overall frequency localization, and an angular one, $\hat{\gamma}(\bar{\xi})$, which controls the directionality, with $\xi$ being the frequency variable and $\bar{\xi}=\xi/|\xi|$. The mother wavelet is thus given by $\hat{\psi}({\xi})=\hat{\gamma}(\bar{\xi})\,\hat{h}(|{\xi}|)$, cf. Fig. 3, left. The whole family of functions used to represent arbitrary signals is then generated by dilation by $2^{j}$, $j\geq 0$, translation by $k\in\mathbb{Z}^{2}$, and rotation by $\theta_{t}$ with $t\in\{0,\cdots,T_{j}\}$.

The angular window $\hat{\gamma}(\bar{\xi})$ is conveniently described using a Fourier series in the polar angle $\theta_{\xi}$, i.e $\hat{\gamma}(\bar{\xi})=\sum_{n}\beta_{j,n}^{t}\,e^{in\theta_{{\xi}}}$. In the frequency domain, a polar wavelet is thus given by

[TABLE]

with the Fourier series coefficients $\beta_{j,n}^{t}$ controlling the angular localization. In the simplest case, $\beta_{n}$ is the Kronecker delta $\delta_{n0}$ and one has isotropic, bump-like wavelet functions, cf. Fig. 2, left. Conversely, when the support of the $\beta_{j,n}^{t}$ is over all integers $\mathbb{Z}$ then one can describe angular windows that are compactly supported in the polar angle $\theta_{\xi}$. Eq. 1 then encompasses ridgelets (Candès and Donoho, 1999; Donoho, 2000) and second generation curvelets (Candès and Donoho, 2005b), cf. Fig. 3, right, which provide quasi optimally sparse representations for image-like signals.

A beneficial property of polar wavelets is that their spatial representation, given by the inverse Fourier transform of Eq. 1, can be computed in closed form. Using the Fourier transform in polar coordinates one obtains (Lessig, 2018b)

[TABLE]

where $h_{n}(|{x}|)$ is the Hankel transform of $\hat{h}(|{\xi}|)$ of order $n$. For $\hat{h}(|{\xi}|)$ we employ the window proposed for the steerable pyramid (Freeman and Adelson, 1991; Portilla and Simoncelli, 2000) since $h_{n}(|{x}|)$ then has a closed form expression (Lessig, 2018b). Note that the angular localization of the wavelets, which is described by the $\smash{\beta_{j,n}^{t}}$ in Eq. 1 and Eq. 2, is invariant under the Fourier transform and only modified by the factor of $i^{n}=e^{in\pi/2}$ that implements a rotation by $\pi/2$.

As indicated in Fig. 3, left, the radial mother window $\hat{h}(|\xi|)$ is defined away from the origin and dilation by $2^{-j}$ centers it at higher and higher frequencies as $j\geq 0$ grows. To represent arbitrary signals, including those whose Fourier transform is nonzero around the origin, an additional window $\hat{g}(|\xi|)$ is required that has support in the disk-neighborhood around $\xi=0$ (and complements $\hat{h}(|\xi|)$). The translates of the inverse Fourier transform of $\hat{g}(|\xi|)$ yield the so called scaling functions $\phi_{k}({x})$ and in our case these will always be isotropic. To simplify notation we will write $\psi_{-1,k}({x})\equiv\phi_{k}({x})$; we refer to (Daubechies, 1992, Ch. 5) for more details on the concept of scaling functions.

The polar wavelets in Eq. 2 together with the just defined scaling functions provide a Parseval tight frame for $L_{2}(\mathbb{R}^{2})$. Thus any function $f({x})\in L_{2}(\mathbb{R}^{2})$ can be represented as (Unser and Chenouard, 2013)

[TABLE]

where $\langle\,,\rangle$ is the $L_{2}$ inner product. Although the above frame is redundant, since it is Parseval tight it still affords most of the conveniences of an orthonormal basis, e.g. the primary and dual frame functions coincide and the norm of the signal equals those of the expansion coefficients. For an isotropic frame the redundancy is thereby $1+1/4+1/4^{2}+\cdots=4/3$ and it increases when directional basis functions are used, i.e. with $T_{j}>1$. As for curvelets, where $T_{j}$ grows according to a parabolic scaling law as a function of $j$, anisotropic representations are typically sparser, which can compensate for the larger redundancy.

Since the definition of polar wavelets in the frequency domain uses a compactly supported radial window $\hat{h}(|\xi|)$, the wavelets have non-compact support in space. Thus, a finite signal representation is not exact, since it requires a truncation of the basis functions. Nonetheless, with a sufficiently large apron region around an image, an arbitrary accuracy can be attained; in our experiments an apron of $4$ pixels sufficed to meet the requirements of photographic applications. To obtain the wavelet representation of a signal, we compute it using a coarse-to-fine, fast wavelet transform-like algorithm where on each level the computations can either be performed using discrete convolutions with filter taps in the spatial domain or by multiplication with the windows in the frequency domain.

To simplify notation, we will typically employ the multi-index $s=(j,k,t)$ introduced in Eq. 1 and, when confusion might arise, write $s=(j_{s},k_{s},t_{s})$. The index $s$ runs over the set $\mathcal{S}$ that a priori includes all scales, translations, and orientations. The cardinality of a set will be denoted by $|\cdot|$, e.g. $|\mathcal{S}|$.

4. A Sheared Local Fourier Slice Equation for Computational Imaging

In this section we derive local Fourier slice photography, our image reconstruction technique from a sparse polar wavelet representation of a light field. We begin by fixing notation and recalling the image formation model. At the end we analyze the computational complexity of the technique as well as sources of error.

4.1. Image Reconstruction Model

We use the two-plane parametrization (Chai et al., 2000) for the light field $\ell(x,y,u,v)$ with $(x,y)$ being the coordinates on the image plane and $(u,v)$ those on the lens, see Fig. 5 for a schematic depiction. We also assume $\ell(x,y,u,v)$ is non-zero only over the camera sensor in $(x,y)$ and over the lens in $(u,v)$ and that it already contains the $\cos{(\theta)}^{4}/F^{2}$ foreshortening factor, where $\theta$ is the angle between the ray and the image plane normal.

For image formation we use the same model as (Ng, 2005). Hence, the image $I(x,y)$ is obtained from the light field as

[TABLE]

where $\alpha=F/F^{\prime}$ is the refocusing parameter, with $F$ and $F^{\prime}$ being the original and new distance of the sensor to the lens plane, respectively, cf. Fig. 5. By changing the notation for the light field to $\ell(x,u;y,v)$, the last equation can be written more compactly as

[TABLE]

where the shear matrix $S_{\alpha}$ is given by

[TABLE]

The shear $S_{\alpha}$ amounts to the transport of the light field in the camera from the original sensor location to the refocused one (Chai et al., 2000; Durand et al., 2005).

To obtain a form of $\ell(x,u;y,v)$ that is amenable to compression, e.g. by thresholding small coefficients, we represent it in the polar wavelets introduced in Sec. 3. Respecting the separable structure of the refocusing in Eq. 4, that is performing one transform over $x$-$u$ and a second one over $y$-$v$, we obtain

[TABLE]

where $\mathcal{L}$ is the index set for the representation that runs over all coefficients $s=(j_{s},k_{s},t_{s})$ and $r=(j_{r},k_{r},t_{r})$. Inserting this representation into Eq. 4 yields

[TABLE]

The last equation shows that it suffices to determine the effect of the sheared projection for the basis functions $\psi_{s}$ and $\psi_{r}$ and that this can be done independently for each of them.

Next, we will thus study general, two-dimensional sheared projection using polar wavelets. This will yield our sheared local Fourier slice equation. We return to imaging in Sec. 4.3.

4.2. A Sheared Local Fourier Slice Equation

Let $f(x,u)$ be a two dimensional signal. We consider the sheared projection

[TABLE]

where $S_{\alpha}$ is given by Eq. 7. When $f(x,u)$ is represented in polar wavelets the equation becomes

[TABLE]

Inspired by Ng’s work (2005), we will seek a numerically practical solution to Eq. 11 in the Fourier domain; a depiction of our approach is shown in Fig. 2. By the Fourier slice theorem, the integral in the last equation can be written as

[TABLE]

where $S_{\alpha}^{-T}(\xi_{x},\xi_{u})^{T}|_{\xi_{u}=0}$ is a linear slice in $\xi_{x}$-$\xi_{u}$ frequency space. By expanding $\hat{\psi}_{s}$ using the definition in Eq. 1, defining $\xi_{0}=(\xi_{x},0)^{T}$, and writing out the inverse Fourier transform we obtain

[TABLE]

that is, $\zeta_{n}^{s,\alpha}(x)$ is the inverse Fourier transform of the $\xi_{x}$-dependent term. Importantly, for our choice of the radial window the $\zeta_{n}^{s,\alpha}(x)$ have a closed form expression, see Appendix A, and the above equation can therefore easily be realized numerically. Furthermore, the original shift $k_{s}$ of $\psi_{s}$ becomes after slicing $P_{x}(2^{-j_{s}}S^{-1}{k_{s}})$, where $P_{x}$ is the projection onto the $x$-axis. This means the shape of the $\zeta_{n}^{s,\alpha}(x)$ remains independent of $k_{s}$. Eq. 13 furthermore shows that the angular localization coefficients $\beta_{j,n}^{t}$ are invariant under the inverse Fourier transform. The sheared projection of an arbitrary polar wavelet $\psi_{x}(x)$ is thus

[TABLE]

and inserting Eq. 14 into Eq. 11 we obtain

[TABLE]

Eq. 15 is our sheared local Fourier slice equation with the $\zeta_{s}^{\alpha}(x)$ being the reconstruction filters that implement projection directly from the wavelet representation. A simple verification of Eq. 15 for a two dimensional Gaussian, for which the ground truth has a closed form solution, is provided in Fig. 7.

The reconstruction kernels $\zeta_{s}^{\alpha}(x)$ are wavelet-like in that they are compactly supported in the frequency domain around $\xi_{x}=2^{j_{s}}$ and well localized in the spatial domain. The former holds since $\hat{\zeta}_{s}^{\alpha}(\xi)$ is a slice of a compactly supported wavelet centered at this frequency, cf. Fig. 6, and the latter since the sliced window has the same decay as $\hat{\psi}_{s}(\xi)$, see Fig. 4. The wavelet-like properties enable a local, sparse reconstruction with a coefficient $f_{s}$ only having a non-negligible effect to $g(x)$ in a small neighborhood around the projected locations $2^{-j_{s}}P_{x}(S^{-1}k_{s})$. Thus, only the wavelet coefficients defined at locations in a sheared tube in the $u$-direction above a location $x^{\prime}$ contribute to $g(x^{\prime})$, see Fig. 9 left. We can therefore think of Eq. 15 as a wavelet representation of the projected signal with the wavelets $\zeta_{s}^{\alpha}(x)$ which are located at non-canonical locations $P_{x}(2^{-j_{s}}S^{-1}k_{s})$. In fact, with isotropic wavelets and $\alpha=1$ the projection in Eq. 15 becomes

[TABLE]

which is a standard, one dimensional wavelet representation of $g(x)$. Our result then also coincides with those of (Lessig, 2018a).

Next to the spatial position of $\psi_{s}(x)$, the contribution of a coefficient $f_{s}$ to the projected signal $g(x)$ also depends on the orientation of $\psi_{s}(x)$, or, equivalently, on the corresponding $\smash{\beta_{j,n}^{t}}$, cf. Eq. 14. As is apparent from Fig. 6, the magnitude of the reconstruction filters $\zeta_{s}^{\alpha}(x)$ is non-negligible only when the direction $S_{\alpha}^{-T}\xi_{0}$ overlaps the effective support of the wavelets in the frequency domain. When curvelet- or ridgelet-like wavelets are used, i.e when $\hat{\psi}_{s}(\xi)$ has strong directional localization, then only a sheared wedge or a small number of wedges from the polar tiling of the frequency plane have effective support over the direction. Hence, only these directions need to be considered in the sum over $s$ in the sheared Fourier slice equation in Eq. 15.

Our result relies on the use of polar wavelets that, through their definition in polar coordinates in frequency space, are compatible with the intrinsic structure of the projection, i.e. with a restriction to a line through the origin in the Fourier domain. With tensor product wavelets, e.g. using Daubechies-type discrete wavelets, one would have a different, skew slicing through the axis-aligned frequency window for every $\alpha$. The reconstruction kernel would then not have a closed form solution and, since the wavelets have no simple description, even determining them numerically would be difficult.

4.3. Imaging with the Sheared Local Fourier Slice Equation

We now return to refocused image reconstruction from the polar wavelet representation of a light field. Using the local Fourier slice equation and inserting Eq. 9 into Eq. 14 we obtain

[TABLE]

The $\zeta_{s}^{\alpha}(x)$ thus provide the reconstruction filters that implement $\alpha$-refocused image reconstruction directly from the wavelet coefficients $\ell_{sr}$ of the light field, that is without the need to obtain a dense pixel representation. In practice, a displayable representation of an image is obtained by evaluating Eq. 17 for every pixel (possibly with multiple samples to increase the quality). We summarize the computations that are required for image reconstruction using Eq. 17 in Algorithm 1.

Eq. 17 can also be thought of as a separable wavelet representation of the image with the wavelets $\zeta_{s}^{\alpha}(x)$ and $\zeta_{r}^{\alpha}(y)$ located at the non-dyadic locations $P_{x}(2^{-j}S^{-1}k)$. Note that the spatial and directional locality $\zeta_{s}^{\alpha}(x)$ discussed in Sec. 4.2 immediately carries over to Eq. 17 and, for example, a coefficient $\ell_{sr}$ contributes to the image $I(x,y)$ only when $\hat{\psi}_{s}(\xi_{x},\xi_{u})$ and $\hat{\psi}_{r}(\xi_{y},\xi_{v})$ are oriented along the slicing direction, cf. again Fig. 6.

4.4. Compressed Image Reconstruction

It is often useful to directly determine the compressed representation of an image from a compressed light field, i.e. without first having to obtain a dense pixel representation as an intermediate step. We will show next how the result of Sec. 4.2 can be extended towards this end.

We assume that the light field is again provided in the polar wavelet representation in Eq. 8 so that the projected light field is given by Eq. 17. Assuming a separable wavelet basis $\psi_{q}^{1}(x)\,\psi_{p}^{1}(y)$ is used for the image, the corresponding expansion coefficients are given by

[TABLE]

By introducing

[TABLE]

we can write this more compactly as

[TABLE]

When $\psi_{q}^{1}(x)$ and $\psi_{p}^{1}(y)$ are one dimensional, bandlimited “polar” wavelets that use the same radial window $\smash{\hat{h}(\xi_{x})}$ as in Sec. 3, then the $\gamma$-coefficients in Eq. 19 can be computed in closed form; the expressions are provided in the accompanying Mathematica code. By the compact support of the wavelets in the frequency domain, the $\gamma_{sq}$ are then, for moderate $\alpha$, non-negligible only when $\max(0,|j_{q}-1|)\leq j_{s}\leq|j_{q}+1|$; when $|\alpha-1|$ is large then $\max(0,|j_{q}-2|)\leq j_{s}\leq|j_{q}+2|$ holds. A plot of the $\gamma_{sq}$ as a function of $P_{x}(2^{-j}S^{-1}k_{s})$ for different values of $\alpha$ and $j_{s}=j_{q}$ is provided in Fig. 8. As can be seen there, the coefficients have a shape similar to $\zeta_{s}^{\alpha}(x)$, since $\hat{\zeta}_{s}^{\alpha}(\xi)$ is essentially a smoothed box function, and in particular they have the same spatial decay.

A consequence of Eq. 20 is that sparsity in the wavelet representation of a reconstructed image is induced by sparsity in those of the light field. In particular, for fixed $(q,p)$ only the $\ell_{sr}$ with $|k_{q}-P_{x}(S_{\alpha}^{-T}k_{s})|\lesssim 2^{-j_{q}}$ and $|k_{p}-P_{x}(S_{\alpha}^{-T}k_{r})|\lesssim 2^{-j_{p}}$ contribute to $\ell_{qp}$, since the wavelets decay in space with dyadic dilation from level to level. Hence, when all $\ell_{sr}$ in the sheared tube above $k_{q}$ and $k_{p}$ are negligible, cf. Fig. 9, then also $\ell_{qp}$ is negligible. The coefficient $\ell_{qp}$ can also become small through cancellation, since both the $\ell_{sr}$ and $\gamma_{sq}$ are signed. Intuitively, this happens, for example, when a region is defocused in the sheared projection and hence the wavelet coefficients on fine levels there have to be small. Such decay estimates are beyond the scope of the present paper and will be investigated elsewhere; existing result in this direction can be found in (Quinto, 1993, 2007).

4.5. Computational Complexity of Image Reconstruction

In the following, we will analyze the computational complexity of Algorithm 1. We assume that the reconstruction is performed for a region $A\subseteq[0,1]^{2}$, which is potentially a subset of the normalized original image plane $[0,1]^{2}$, with size $|A|$ and that $R$ is the sampling rate per pixel in the reconstructed image. The input light field is assumed to have resolution $N_{x}\times N_{x}\times N_{u}\times N_{u}\times 3$ (we will ignore the apron that is used in practice to avoid boundary artifacts, since it is small and has a negligible effect on the complexity yet would make the analysis considerably more cumbersome).

Step 1 and Step 2 in Algorithm 1 have negligible costs and will be disregarded. For Step 3, the separable reconstruction kernel $\zeta_{s}^{\alpha}(x)\,\zeta_{r}^{\alpha}(y)$ has effective support $2^{-j_{s}}W\times 2^{-j_{r}}W$, where $W$ is a window function-dependent constant (see Fig. 4). The computational costs $K_{sr}$ that arise for each coefficient $\ell_{sr}$ are thus

[TABLE]

where $c_{\zeta}$ is the cost for evaluating $\zeta_{s}^{\alpha}(x)$ at one point. The number of coefficients on level $(j_{s},j_{r})$ in a dense wavelet representation is

[TABLE]

since at each level the resolution is reduced by a factor of $2$ in each direction. As before, $T_{j}$ is the number of different orientations on level $j$ and the factor of $4=3\cdot 4/3$ accounts for the three color channels and the redundancy of the isotropic frame. The computational costs K_{j_{s},j_{r}}=K_{sr}\,\big{|}\mathcal{L}_{j_{s},j_{r}}\big{|} for image reconstruction on level $(j_{s},j_{r})$ are for a dense polar wavelet representation thus

[TABLE]

For the costs $K_{j_{max}}$ across all levels we then have

[TABLE]

where $P_{j_{max}}=2^{-4{j_{max}}}(2^{{j_{max}}+1}\!-\!1)^{2}$.

Two facts can be exploited to reduce the costs $K_{j_{max}}$. First, with anisotropic, curvelet-like frame functions only the $\alpha$-sheared ones overlapping the slicing direction are required, cf. Fig. 6. These are of order one so that the costs $K_{{j_{max}}}^{\alpha}$ are

[TABLE]

where $c_{T}=\mathcal{O}(1)$. The increased redundancy of a directional representation system hence does not affect the costs since the number of frame functions overlapping the slicing direction remains constant. The second reduction of the costs $K_{j_{max}}$ arises from the sparsity in a light field’s wavelet representation, i.e. that only a sparse set $\smash{\mathcal{L}_{j_{s},j_{r}}^{\alpha,\epsilon}}$, consisting of non-negligible ones with respect to a compression parameter $\epsilon$ (in the simplest case of hard thresholding, $\epsilon$ gives the threshold), suffices to represent the signal. The ratio between the number of coefficients in the full and sparse representations is known as compression rate

[TABLE]

Assuming it is independent of the level, i.e. $\mathrm{cr}_{j_{s},j_{r}}^{\epsilon}=\mathrm{cr}_{\epsilon}$, when sparsity and directionality are exploited the costs $K_{{j_{max}}}^{\alpha,\epsilon}$ then are

[TABLE]

With basis dependent constants being ignored, this becomes in big-O notation

[TABLE]

This shows that the costs scale linearly in the area that is to be reconstructed and the sampling rate and as $1/cr^{\epsilon}$ in the compression rate. Thus, as the number of coefficients in the sparse representation $\mathcal{L}_{j_{s},j_{r}}^{\alpha,\epsilon}$ decreases and the compression rate grows also the computational costs decrease.

The complexity of Ng’s Fourier slice photography (Ng, 2005) is, in our notation, \mathcal{O}\big{(}|A|\,R\,N_{x}^{2}\big{)}. Although both works exploit that the projection becomes trivial in the Fourier domain, the spatial localization of the wavelets in our approach introduces again the directional resolution parameter $N_{u}$ in Eq. 28. However, with the wavelets we also have a dependence on the compression rate $\mathrm{cr}^{\epsilon}$. This can compensate for the $N_{u}$ factor, depending on the rate $\mathrm{cr}^{\epsilon}$ that can be attained for the light field. Although theoretical characterizations of $\mathrm{cr}^{\epsilon}$ exist, see e.g. (Mallat, 2009, Ch. 6, Ch. 9), these require technical assumptions about the signal that are difficult to precisely meet in practice. We will hence not pursue a further theoretical analysis here. Nonetheless, the practical utility of wavelets for the compression of image like signals, and hence that significant compression rates can be attained, is by now well established, as is evidenced by their use in the JPEG2000 standard.

4.6. Image Reconstruction Error

Ng (2005) discusses two sources of error for image reconstruction, namely roll off error and aliasing. We avoid roll-off error since our reconstruction kernels $\zeta_{s}^{\alpha}$ are the exact ones for our analysis wavelets, which in turn provide a Parseval tight frame; when the wavelets are used up to level $j_{\mathrm{max}}$ then all signals with bandlimit $2^{j}\pi$ can be represented exactly. Aliasing is of concern because the shearing can increase the bandlimit (the change in the bandlimit can be seen, e.g. in bottom right plot in Fig. 6). In our approach, this implies that for $\alpha<1$ the $\zeta_{s}^{\alpha}$ have a bandlimit beyond those of the original polar wavelets $\psi_{s}$. Consequently, they also need to be defined over a finer grid than the original wavelets to allow for perfect reconstruction. But, as shown in Fig. 9, left, the shearing also affects the grid over which the basis functions are defined and through this the $\zeta_{s}^{\alpha}$ are inherently defined on a lattice that has the appropriate density.

Beyond roll off error and aliasing, a third source of error in fact arises in Ng’s work (2005). As depicted in Fig. 9, right, the discrete Fourier transform yields a lattice of discrete frequencies in the Fourier domain. Except in the trivial case when $\alpha=1$, the slicing direction $S^{-T}\xi_{0}$ does not lie on the lattice and the resampling onto a regular grid along the slicing direction provides an additional source of error. This error does not occur in our approach since our basis functions are defined in the continuous Fourier domain and we analytically compute the projection as a function of $\alpha\in\mathbb{R}$.

The main source of inaccuracies in our approach are in practice those introduced by a finite truncation of the basis functions or, equivalently, of the filter taps used in the fast transform. As discussed in Sec. 3, these can be ameliorated by using appropriate padding and filter sizes.

5. Experiments

In this section we present experimental results that verify the practicality of the image reconstruction technique developed in the last section. Additional results as well as the raw images for most of the presented figures are provided in the supplementary material.

5.1. Data sets

We implemented a custom light field film class for the pbrt renderer (Pharr and Humphreys, 2010) so that we could easily vary the spatial and angular resolutions as well as the optical properties of light field data sets. The film class directly records light fields in the two plane parameterization, applies the foreshortening factor of $\cos{\theta}^{4}/F^{2}$ and also computes a depth map for the scene. We used pbrt to generated synthetic light fields for a dragon scene and a villa interior, both with a resolution of $1025\times 1025\times 33\times 33\times 3$. We also performed experiments with photographic light fields from the Stanford Lytro light field archive.111http://lightfields.stanford.edu/LF2016.html We selected 20 light fields from different categories, all with a resolution of $541\times 376\times 14\times 14\times 3$.

All data sets were high dynamic range with computations performed in single precision. We used two to three wavelet levels in all experiments (because of the limited resolution in the angular dimension) and for the experiments with the photographic light fields the same parameters (such as range of $\alpha$ values and tone mapping parameters) were applied in each case, which might be sub-optimal in individual instances. Brightness variations that can be seen in some of the results for varying $\alpha$ stem from (independent) tone mapping.

5.2. Implementation

We developed a Mathematica reference implementation of Algorithm 1 (available in the supplementary material) and a basic, multi-threaded C++ implementation. Since the reconstruction filters $\zeta_{s}^{\alpha}(x)$ are relatively expensive to evaluate, cf. Appendix A, we sampled them in a preprocessing step and interpolated at runtime (the error introduced through the interpolation was below $10^{-7}$ and thus negligible for photographic applications). Post-processing bias was avoided by using one sample per pixel for image reconstruction and no interpolation filter on the image plane. “Reference” solutions were similarly computed using naïve projection in the pixel domain without filtering of the light field data sets or the pixel data.

For the polar wavelets, we used filter taps of size $81\times 81$ and, as mentioned earlier, an apron of 4 pixels. Larger values did not improve the reconstruction. Note that our algorithm for refocused image reconstruction is itself parameter free.

To study the effect of sparsification on image reconstruction quality we implemented a simple transform coding scheme with hard thresholding, i.e. we set to zero all coefficients whose magnitude is below a threshold $\epsilon$ dependent on the light field’s $L_{2}$ norm to compensate for overall brightness differences between the data sets. Results will be reported using either the number of nonzero coefficients in the light fields $\epsilon$-sparse wavelet representation (denoted by nzs) or the compression rate (denoted by cr), i.e. the number of nonzero coefficients in a sparse over the total number in a dense wavelet representation. The number of nonzero coefficients provides an indication of the storage requirements although it is an upper bound since our wavelet representation only provides the transform coding step of a full compression scheme and we do not consider blocking, quantization, entropy coding and other aspects that are critical in practical compression algorithms.

5.3. Simple experiments

To demonstrate the correctness of the sheared local Fourier slice equation as well as to gain some understanding of various conceptual aspects we performed experiments on two dimensional signals yielding a one dimensional projection.

Basic verification

We verified the correctness of the local Fourier slice equation using the sheared projection of a two dimensional Gaussian for which an analytic solution exists. As shown in Fig 7, our reconstruction matches the analytic one very well with a maximum error on $1.86\times 10^{-6}$. This is of the same order as the reconstruction error of the 2D input signal, and hence attributable to inaccuracies in the transform yielding the wavelet coefficients.

Non-smooth signals and effect of $\alpha$

To obtain insights on the behavior of our technique for “natural images” as well as to understand the effect of $\alpha$ on the reconstruction quality we considered the projection of a monochromatic dragons scene image (derived from the image in Fig. 12). The results in Fig. 11, right, demonstrate that the error is smaller than what can be perceived visually and on the same order as differences resulting with different reconstruction kernels for the pixel domain projection (not shown). Shown in the figure is also a quantitative analysis of the error as a function of $\alpha$, demonstrating that only a mild dependence on the angle exists. The results, furthermore, reveal that our technique provides slightly lower error rates than classical Fourier slice photography (Ng, 2005). The plot on the top left of Fig. 11 indicates that the higher errors for the technique result from the interpolation from the axis-aligned DFT grid that is required for it, see Fig. 9. Since our wavelets are defined in the continuous Fourier domain and the restriction to a slice is performed there, no such interpolation is required.

Sparsity and Directionality

For the dragon scene image we also studied the effect of sparsity as well as the gains that are possible using oriented, curvelet-like frame functions. The left plot in Fig. 10 shows the error as a function of the compression rate with hard thresholding. As one would expect for a wavelet representation, very accurate reconstructions are possible with a small fraction of the full coefficient set. Furthermore, the error increases smoothly with the compression rate. The middle plot in Fig. 10 depicts the norm $\|\zeta_{jkt}(x)\|$ of the reconstruction kernel $\zeta_{jkt}(x)$ as a function of the orientation $t$ of the wavelets. Because of the localized support of the directional wavelets in the frequency domain, cf. Fig. 6, only $\zeta_{jkt}(x)$ whose orientation $t$ matches the slicing direction $\smash{S_{\alpha}^{-T}\xi_{0}}$ are non-negligible. The right plot in Fig. 10 verifies that orientations far from the projection direction do not have to be taken into account for reconstruction, i.e. the reconstruction error is sufficiently small when these are ignored. Furthermore, since for typical values $\alpha$ does not fundamentally change the direction, some orientations are irrelevant, e.g. $t\leq\pi/2$ and $t\geq 3\pi/2$ in the plot. Thus, the frame coefficients for the negligible orientations do not have to be stored, even if these are above an $\epsilon$-threshold. $\|\zeta_{jkt}(x)\|$ can hence be understood as a signal independent form of sparsification.

5.4. Image reconstruction

Refocusing

The reconstruction of a refocused image is provided in Fig. 12. The comparison to the projection in the pixel domain, also shown in the image, verifies that our result is visually indistinguishable and that the differences between both are on the order of what one would obtain with different reconstruction kernels for projection in the pixel domain.

Fig. 13 shows reconstructions with only a subset of levels and Fig. 19 the contribution made by individual levels. It can be seen that the error increases gracefully as one disregards finer levels with the reconstructed images becoming smoother but largely without objectionable artifacts, though some ringing can be observed for $j\leq-1$ and $j\leq 0$ in the in-focus region. Since a reconstruction with a subset of levels is considerably cheaper, requiring only $0.098\%$, $1.72\%$, $13.43\%$ of the computation time for all levels for $j\leq-1$, $j\leq 0$, and $j\leq 1$, respectively, substantial savings for out-of-focus regions are possible when a depth map is known.

Sparsity

Fig. 14 and Fig. 15 show the relationship between the number of nonzero coefficients and $L_{1}$ reconstruction errors for 10 different photographic light fields for $\alpha=0.6,0.85,1.05,1.35$ (top to bottom in each figure). For $\alpha=0.6$, in both figures one data set (the fifth in Fig. 14 and the third in Fig. 15) yields considerably larger errors than the other ones. For the two data sets large areas with high frequency foliage, requiring many small wavelet coefficients for an accurate representation, are in focus for $\alpha=0.6$ and this behavior is thus to be expected. A similar observation holds true for the fourth and fifth data set in Fig. 15 although it is less pronounced there since the in-focus regions are smaller. In contrast, for the third data set in Fig. 15 one has, except for some edges, a low visual complexity independent of $\alpha$. Correspondingly, one has a consistently low reconstruction error. For the remaining data sets there is for each $\alpha$ one region with high visual complexity in focus and the reconstruction errors remain thus relatively constant. Qualitatively equivalent results hold for the $L_{2}$ and $L_{\infty}$ norms. Results for 10 different photographic light fields can be found in the supplementary material.

Fig. 18 provides a visual comparison of reconstructed images as the sparsity in the light field’s wavelet representation increases. The figure demonstrates that the image fidelity degrades gracefully as the compression rate increases. Furthermore, when the error becomes visible then it amounts to a lack of high frequency features but largely without objectionable artifacts (image sequences as a function of the compression rate are provided in the supplementary material). The storage requirements for the coefficients of the compressed light fields in Fig. 18 (with $cr>1$) are $2.23$ GB, $1.27$ GB, $660$ MB, $122.39$ MB, and $40.23$ MB, respectively. As a comparison, the dense light field, required for projection in the pixel domain or Fourier slice photography, requires 13 GB of storage and the uncompressed wavelet representation 48 GB.

In Fig. 16 we report the $L_{2}$ error in the reconstructed image as a function of the number of nonzero coefficients for our local Fourier slice photography and a JPEG-compressed representation of the light field (similar to (Alves et al., 2018)). For the latter, the sparse representation was obtained by considering each slice in the light field (for fixed u-v index) as an image and applying the transform coding step of JPEG, consisting of mask-weighted quantization in the discrete cosine transform (DCT) domain over $8\times 8$ blocks. We did not apply JPEG’s entropy coding step to obtain a roughly fair comparison to our technique that also only implements transform coding. Fig. 16 demonstrates that for large coefficient counts, i.e. little compression, classical JPEG outperforms our approach. This is to be expected since polar wavelets are redundant and hence more data is required for perfect reconstruction when compared to the non-redundant DCT used in JPEG. For low coefficient counts, i.e. more aggressive compression, local Fourier slice photography provides, however, lower errors than the JPEG-compressed representations. Furthermore, and this is at the heart of our work, for the JPEG-compressed light field data sets a dense representation has to be obtained for refocused image reconstruction, requiring GBs of memory, while our approach can perform the projection directly from the sparse wavelet representation.

All-in-focus reconstruction

Fig. 20, right, shows a reconstruction of the dragons scene as well as of the checkerboard ground plane with a depth dependent $\alpha$ value so that the entire scene is in focus. Slight artifacts are visible around the dragon silhouettes, since we do not take the varying support of the $\zeta_{s}^{\alpha}(x)$ into account and only sample the depth map at the location of the reconstruction kernels.

Performance

Fig. 17 shows the relative execution time as a function of the compression rate. Although our implementation is not particularly optimized and we only use the Eigen library for the sparse wavelet representation, the results demonstrate that sparsity can lead to a substantial reduction in computation time. The plots also show the execution time decreases approximately as $1/\mathrm{cr}$, as one would expect from our analysis of the computational complexity in Sec. 4.5, see in particular Eq. 28.

The absolute computation time of image reconstruction is currently approximately two minute for a $1025\times 1025\times 33\times 33\times 3$ light field data set with a dense wavelet representation and $13$ seconds for the highest compression rate we considered. The projection in the pixel domain requires in our implementation $25$ seconds from a decompressed, dense light field and roughly the same time is required for the sparse projection with a compression rate of $200$. The computation of the wavelet representation of the light field requires approximately a minute.

5.5. Discussion

Our experimental results demonstrate the practical viability of Algorithm 1 for the reconstruction of refocused images from the sparse wavelet representation of a light field. We verified that high fidelity images can be obtained from a highly sparse representation and that the error increases gracefully with the compression rate. Furthermore, our experiments show that the error depends on the visual complexity of the in-focus region, which can be exploited when a depth map is available. Additionally, we demonstrated that simple hard thresholding of polar wavelet coefficients is, at least for moderate to high compression rates, competitive with the transform coding step of JPEG, which uses highly optimized masks.

Our results also show that sparsity in the wavelet representation can lead to a reduction in the computation time through the smaller number of coefficients that has to be processed, although our implementation is currently slower than projection in the pixel domain when a dense representation of a light field is directly available. Algorithm 1 is easily parallelized by exploiting that the reconstruction for each pixel is independent, i.e. there is $N_{x}\times N_{y}$ data parallel work. This makes it directly amenable to a GPU implementation where one could also take advantage of half-precision, which is sufficient to obtain artifact free images. This would also provide benefits on the embedded processors typically available in cameras.

Fig. 20 showed first results on the reconstruction of all-focus images from the compressed wavelet representation. While currently not artifact free when the depth map contains discontinuities, the results verify the potential of our approach to obtain all-focus images, which is not possible using Fourier slice photography where a fixed $\alpha$ has to be used. To remove the current artifacts, the depth map needs to be preprocessed in a mip-map-like manner so that an average depth can be sampled at each wavelet level and it might also be necessary to restrict the support of wavelets at depth discontinuities.

6. Future Work

The local Fourier slice photography algorithm developed in the foregoing suggests many directions for future work.

In our work, we consider the sparse wavelet representation of light fields, which corresponds to the transform coding step of a compression technique. To make local Fourier slice photography practical for applications, this has to be extended to a full compression scheme by also considering, for example, blocking, quantization, and color coding. Our comparison to the transform coding step of JPEG indicates already that the sparse representation will translate to significantly reduced storage requirements also in practice.

Our technique would also benefit from additional work on the wavelet representation. For example, the radial window $\hat{h}(|\xi|)$ we currently employ does not provide very good decay in the spatial domain and, similar to (Ward et al., 2015), one could investigate how better radial windows can be constructed. One should also further investigate the trade off between increased sparsity for directional, curvelet-like representations and the higher computational costs for evaluating these. Our results in Sec. 4.4 showed that sparsity in the light field induces sparsity in the reconstructed image. A quantitative description of this could potentially lead to a further reduction of the storage requirements as well as computation times.

By representing the light field in polar wavelets, many existing techniques for editing and processing are no longer directly available. However, we believe that many of them can be translated to the polar wavelet domain and it might, in fact, provides advantages to conventional approaches. For instance, the feature-aware resizing of Gastal and Oliveira (2017) is naturally formulated in polar coordinates and could hence be performed directly in a sparse representation. We believe that other tasks, such as the shearlet-based light field reconstruction from a sparse set of views in (Vagharshakyan et al., 2018) or learning based techniques such as (Kalantari et al., 2016) could also benefit from the efficacy of a polar wavelet representation.

In the present work we considered refocused image reconstruction. Another important application of light fields is novel view synthesis. It would be interesting to investigate if this can also be performed directly from a sparse polar wavelet representation.

Our current approach for image reconstruction exploits the separability of the refocusing problem so that two dimensional polar wavelets are sufficient. This is sub-optimal concerning the compression rates that can be attained for the light field and also since one obtains a separable wavelet representation for the reconstructed image. With 4D polar wavelets, which can be constructed as an extension of the polar wavelets used in the present work, cf. (Ward and Unser, 2014), the entire light field could be represented in one wavelet basis and the projection of the data would then yield two dimensional, curvelet-like polar wavelets. In our opinion, this is both theoretically and practically an interesting direction for future work.

7. Conclusion

In this paper, we presented local Fourier slice photography, an algorithm to reconstruct refocused images from a sparse wavelet representation of a light field. For this, we derived a sheared local Fourier slice equation, which extends the local Fourier slice equation of (Lessig, 2018a) to the case of sheared projection. The equation provides analytic, wavelet-like reconstruction kernels for obtaining a refocused image directly from a light field’s wavelet coefficients. The direct reconstruction from a sparse representation is at the heart of our work and it avoids the need for a dense light field representation that exists for techniques in the literature. We experimentally verified that high fidelity images can be reconstructed from a highly sparse representations of a light field, providing the potential for significant storage requirements, and our results demonstrate that image quality degrades gracefully as the compression rate increases. Furthermore, also the computational costs can be reduced by directly reconstructing images from the sparse wavelet representation. We analyzed this theoretically and demonstrated efficiency gains as a function of compression rate experimentally. Because of the spatial localization of the wavelets, the costs of our technique depend on the light field’s angular resolution. This was not the case for Fourier Slice Photography (Ng, 2005) but the localization enables us, for example, to obtain all-focus images, which is not possible using this technique.

Acknowledgements.

The comments by the anonymous reviewers helped to considerably improve the manuscript, in particular their insistence to try the algorithm on photographic light fields. The Stanford light field group is acknowledged for making the Lytro data sets available to the community. First ideas for the project were developed while the author was a post-doc in Marc Alexa’s computer graphics group at TU Berlin. Many thanks also to Eugene Fiume for continuing support.

Appendix A Spatial representation of reconstruction filters $\zeta_{s}^{\alpha}$

When $\hat{h}(|\xi|)$ is the wavelet window from the steerable pyramid, the $\zeta_{s}^{\alpha,m}(x)$ in Eq. 14 are given by

[TABLE]

where

[TABLE]

and

[TABLE]