Efficient computation of backprojection arrays for 3D light field deconvolution

Light field deconvolution allows three-dimensional investigations from a single snapshot recording of a plenoptic camera. It is based on a linear image formation model, and iterative volume reconstruction requires to define the backprojection of individual image pixels into object space. This is effectively a reversal of the point spread function (PSF), and backprojection arrays H' can be derived from the shift-variant PSFs H of the optical system, which is a very time consuming step for high resolution cameras. This paper illustrates the common structure of backprojection arrays and the significance of their efficient computation. A new algorithm is presented to determine H' from H, which is based on the distinct relation of the elements' positions within the two multi-dimensional arrays. It permits a pure array re-arrangement, and while results are identical to those from published codes, computation times are drastically reduced. This is shown by benchmarking the new method using various sample PSF arrays against existing algorithms. The paper is complemented by practical hints for the experimental acquisition of light field PSFs in a photographic setup.


Introduction
Following the seminal publication by L et al. [1], a considerable amount of work has been dedicated to the development of 3D diagnostics for microscopes, using digital plenoptic cameras as imaging devices. Compared to standard photographic cameras, these instruments allow to measure not only the intensity distribution at their image plane, but also record additional directional information on the light rays in the scene. This is realized by inserting an array of microlenses (MLA) into the optical path, close to the image sensor, which distributes light according to its direction onto different pixels. The captured image is therefore a coded representation of the 4D light field, and its spatio-angular information allows to derive depth coordinates, and, for transparent objects, the 3D intensity distribution within the volume. This is a scanless technique based on a single snapshot recording, with the attractive potential of investigating dynamic processes in three dimensions. In its original implementation, light field microscopy mimicked the work flow of traditional 3D deconvolution: Here the input data is in the form of a focal stack, recorded as a sequence while sweeping the object along the optical axis [2]. Image formation is modeled as a convolution of the object space intensity distribution with the point spread function (PSF) of the microscope, which defines light transport within the optical system. Deconvolution methods seek to revert this process, using a known PSF as a tool for estimating the original volume from recorded image stacks. With captured light field data, however, focal stacks can be computed from a single exposure by synthetic refocusing [3], eliminating the time-consuming acquisition of sequences. Capturing of transient processes is therefore only limited by the frame rate of the camera. This comes at a cost: As a plenoptic camera typically has to sacrifice lateral resolution for the additional directional information, early light field microscopy suffered from comparably low pixel counts of the reconstructed volumes. A significant improvement was published by Broxton et al. [4], which circumvented the generation of focal stacks, but instead directly operated on light field data recorded by a plenoptic camera. So-called light field deconvolution is closely related to super-resolution approaches [5], and models the formation of the camera's raw images, i.e. the captured 4D light field, by applying a measurement matrix H to the intensity distribution within object space. For a discretized volume, this matrix, also termed light field PSF, defines the intensity received by each single pixel from all individual light-emitting voxels. These interrelations are sketched in Fig. 1, which also serves as a reference for the definitions and notations that are going to be used throughout the paper. Similar to a standard PSF, H is found by measuring or simulating the sensor's pixel response to a subresolution, point-like light source at all voxel positions. The sketch shows light emitting voxels at two distinct positions within the object space g, and the resulting pixel intensities within the image space f. This technique provided a considerable boost in volume resolution and allowed the 3D analysis of live animals under a microscope [6]. However, resolution was non-uniform across the depth of field and suffered from artifacts especially close to the native object plane. This could be ameliorated by using an enhanced optical setup with additional phase plates, as suggested by C et al. [7]. A different route was taken by S et al. [8], who incorporated a depth-dependent filtering in the light field deconvolution algorithm, which resulted in an artifact-free volume reconstruction. A significant computational speed-up and reduction of artifacts was achieved by L et al. [9] by transferring both the PSF H and the measured light field data into phase space and separating the spatial frequencies in the deconvolution step. At the core of these techniques are deconvolution algorithms, which have been in use before, performing de-blurring of two-dimensional images as well as 3D reconstruction from traditional microscopic focal stacks. All of the above mentioned publications rely on the classical Richardson- Lucy scheme, but a variety of other algorithms may be suitable as well, with an overview given e.g. by S et al. [10]. In addition to the PSF H, most of these algorithms require to define light transport through the optical system in the reversed direction: While H describes a forward projection of the volume onto the image plane, we need to formulate a backprojection H of the pixels into object space. It is, however, not trivial to derive this backprojection from a given light field PSF, because the latter is shift-variant, which means that point sources at different voxel positions create different sensor responses, which are commonly stored in a multi-dimensional array H. While the literature gives detailed explanations of the wave optical modeling of the light field PSF, it does not provide recipes for the computation of H from H, which requires an elaborate procedure. The publications [6,8,9] include computer code that also performs this step, but without specific comments on its functionality. Moreover, computation of H using these codes is slow for modern plenoptic cameras with a high pixel count and MLAs in a hexagonal grid. The contribution of this paper is the following: It first clarifies the significance and structure of H and its relation to the original light field PSF, and complements this by hints on how to record the PSF experimentally in a suitable array. The paper illustrates the relevance of a quick computation of H , and then presents a simple and rapid algorithm for this task. It is purely based on the distinct relation of the elements' positions within H and H , respectively, and a benchmark using sample PSFs shows a substantial speed-up compared to previously published codes while providing identical results.

Image formation and recording of the light field
In a standard photographic camera, single sensor pixels (or grains of a chemical film) integrate the incident light over a certain solid angle. As a consequence, directional information is lost and the images are flat. A plenoptic camera uses a microlens array (MLA) close to the sensor, as sketched in Fig. 2, which distributes light rays according to their direction onto different sensor pixels [11]. The MLA has the effect of a multiplexer [12], coding the lost directional information into the captured raw image, which is a recording of the 4D light field of the scene. Two points in object space, marked A and B in the figure, generate distinct spot patterns on the image sensor. These patterns are representations of the light field, often termed plenoptic function, which defines the transport of light energy along rays in space. In its simplest form this is a 4D function,  with two spatial coordinates describing position, and two angular coordinates defining orientation. Consequently, the raw image recorded by a plenoptic camera can be interpreted as an overlay of two two-dimensional coordinate systems, sketched on the left of Fig. 3. The gray circles illustrate circular microimages formed behind the lenslets of the MLA, and their position is expressed by an outer coordinate system. Little boxes represent sensor pixels, and their position within the respective microimages is given by an inner system. The exact physical meaning of these systems depends on the type of plenoptic camera: In its initial implementation, sampling of spatial and angular data is strictly separated, with the outer and inner coordinate system defining position and ray orientation, respectively [13]. Angular resolution of the sampled light field is defined by the number of pixels within a microimage, and lateral resolution of the images is given by the number of microlenses. This illustrates the already mentioned trade-off of a plenoptic camera, sacrificing pixel resolution for additional directional information. For a so-called focused plenoptic camera, with a different separation between MLA and sensor, sampling of spatial and directional data is intertwined [14]. This allows an increased lateral resolution at the cost of a more complex processing of the raw data. Regardless of the type of camera, a sensor recording of a single light point, as outlined in Fig. 2, represents a point spread function, a spatio-angular light field PSF, that defines the optical system. Such a PSF may look like in the sketch on the left of Fig. 3. The right hand side of the figure shows actual PSF patterns captured experimentally with a plenoptic camera for points at two different depth positions. Here a so-called multi-focus camera was used, which features an MLA with three different lens types in a hexagonal arrangement [15,16]. The pixel patterns formed by discrete light points are, however, not only a function of the point's depth coordinate, but also depend on its lateral position. This holds true even if the plenoptic camera is used at a microscope with a telecentric objective lens. A single PSF is therefore not sufficient to characterize light transport within the system, but a suitable set of PSF patterns must be used for light field deconvolution.

PSF structure
The coordinate system used in this paper and various dimensions and indices are included in the sketch of Fig. 1. The image f is formed on an , -plane and is discretized into × pixels. Note that here the superposition of two coordinate systems is disregarded, but instead absolute pixel positions are used. The volume g is contained within an , , -cube, with the z-axis aligned to the camera's optical axis. It is discretized accordingly into × × voxels. Lateral voxel sizes are tied to the physical pixel pitch of the sensor by the main lens magnification, while voxel depth is arbitrary and subject to the chosen discretization. Light emitted by a single point (or voxel) in space is transferred through object space and the optical system and is then recorded by the sensor pixels as an intensity distribution, as shown in Fig. 1 in blue and red. A complete light field PSF H is therefore a collection of these patterns, captured with a light point at all relevant (voxel-) positions within object space. In its full extend, H establishes a relation between light emitted at each voxel and the received intensity at each pixel. This is a huge amount of data which can hardly be handled by a processing algorithm. Two findings help to drastically cut required information [4]: Point spread patterns of single points are very sparse, so that zero pixels can be largely discarded. This means cropping recorded patterns using rectangular cutouts, illustrated by colored boxes in Fig. 1. And second, the regular arrangement of the lenslets within the MLA results in periodically repeating pixel patterns when shifting a light point in lateral directions. For each axial depth, the optical system is therefore defined by a limited number of light field patterns, which need to be determined with a point source at all voxel positions within a representative area, called an elementary cell. Shape and size of this area depend on the MLA layout, as illustrated in Fig. 4. For lenslets arranged in a rectangular grid, on the left of the figure, the elementary cell is simply quadratic, covering the extends of a microimage. The light field PSF H has to hold sensor responses to a point source placed at all corresponding voxel positions, which is 5×5 in this example, the number of pixels within a microimage. Absolute physical voxel sizes in object space depend on the magnification of the main lens, subject to optional super-or undersampling factors in the code implementation. For a hexagonal MLA layout, shown in the middle of Fig. 4, the elementary cell takes the shape of the dashed box, and shifting the point source beyond the corresponding limits in object space will result in repeating patterns. For the case of a multi-focus camera with three different microlenses, sketched on the right, the elementary cell has to be expanded. It is the smallest possible area of the MLA layout that allows to tile the entire plane with periodic copies. Considering lines of symmetry can help to further reduce its size. This approach is based on the assumption of a shift-invariant PSF of the isolated main lens, which holds, strictly speaking, only for ideal telecentric systems like microscopes [1]. If a plenoptic camera is used in a photographic setup, this assumption may be at stake and must be carefully examined. The general concept of deriving a backprojection array from a given light field PSF is, however, not limited to orthographic projections in microscope systems.
With this background it is clear that image formation on the sensor plane of a plenoptic camera cannot be modeled as a convolution, because the forward projection, defined by H, varies from point to point. In a discretized form, it is given by the generalized linear equation [4] f = H g (1) Image formation can be interpreted as distributing copies of the individual light field patterns of H, scaled by the magnitude of the respective object space voxels in g, which sum up to form the light field image f. In computer code, this distribution can be efficiently implemented as single, discrete convolutions. To be able to do so, it has to be ensured that the patterns are mapped to the correct pixel positions: During acquisition of H, the bounds of the cutouts, drawn as colored boxes in Fig. 1, must be shifted likewise when the point source traverses from voxel to voxel. The light field PSF H has to be given in a structure that is suitable for actual computer applications. Codes published in [6,8,9] conveniently define H as a 5-dimensional array, with two dimensions (s,t) holding the image pixels, and three dimensions (x,y,z) defining the respective object space coordinates. It is therefore a collection of two-dimensional patterns, arranged in a three-dimensional array.

PSF acquisition
For microscope applications, determination of the light field PSF H has commonly been done on a theoretical basis [4,[6][7][8]. This is reasonable, because the optical system is very precise, and due to the high magnification, the required shifts of the point source are on an extremely small scale. For photographic setups, however, it is beneficial to cast this into an experimental calibration procedure to account for imperfections and not precisely defined elements. This requires a very small, point-like light source, close to an isotropic emitter, that is precisely positioned within the calibration volume by a computerized 3-axes translation stage [17]. The stage is synchronized with the plenoptic camera, which records light field raw data at each step, and the cutout pixel patterns are then forged into the array H. Here the cone-like field of view of the photographic main lens needs to be considered, with the consequence that the physical boundaries of the Fig. 5. Summation of all elements of the PSF matrix H within one depth plane. The PSF was recorded experimentally using a photographic main lens, and the structure of the aperture constructed from 9 movable blades is clearly seen.
elementary cells vary with the depth coordinate. It is interesting to examine the summation of all the elements of H within one depth plane, shown in Fig. 5. Here the PSF was acquired experimentally with a photographic plenoptic camera (Raytrix R29 with Nikkor 200mm f/4 main lens) [18]. By subsequently shifting point source and cutout region, the aperture of the main lens is gradually illuminated from all possible directions, filling up the spatio-angular data in the image. This reveals the complex shape of the aperture, designed with 9 movable blades, which are variable and not precisely know, showcasing the requirement of an experimental calibration.

Deconvolution and backprojection array
Volume reconstruction is an inverse problem that seeks to revert Eq. (1), trying to find the original volume g from a measured light field f, using a known PSF H. This is an ill-posed task, and inevitable noise prevents a simple inversion of the image formation equation, so that iterative deconvolution techniques have to be applied instead. One representative of the wide range of different algorithms is the classical Richardson-Lucy scheme, which is still common in advanced deconvolution approaches [19]. It has so far been used in the framework of light field deconvolution, and its iterative update in matrix-vector notation reads [8] Note that this must not be interpreted as standard matrix multiplications, as the PSF is in the form of a multi-dimensional array, which has to be implemented accordingly in the computer code. In essence, the scheme of Eq. (2) computes an error quotient by comparing the measured image f to the forward projection of the current volume estimate, H g ( ) . This error is then backprojected into object space by means of the array H , and updates the volumetric intensity distribution g. This paper now seeks to establish a relation between the light field PSF H and the backprojection array H . The latter has been introduced as a linear operator, which defines a backprojection of a single pixel through the optical system into object space. Alternatively, it shows for every pixel the (voxel-) position and proportion of light in the volume that contributes to the total received intensity [4]. If we think of a volume with unity intensity and consider the common 5-dimensional array form of H, then the forward projection onto the image plane is simply a shifted summation of the individual PSF patterns. This is illustrated in Fig. 6 for a single depth plane and an array H, which holds 3×3 pixel patterns with 3×3 pixels each. Each pattern defines the PSF of a single point at position ( , ) in object space, and is a distribution of pixel intensities in ( , ) coordinates, shown on the left of the figure. Certain elements in the sketch have been given a number and are colored in red or blue. Forward projection of the volume is conceptualized in the middle of the figure. Each of the boxes here represents an image pixel, which receives contributions from several points in space, shown as an overlay of patterns which have been shifted relative to each other according to the point's position. The center pixel intensity, e.g., is a summation of all those labeled elements of H that have been marked with red boxes. Grouping these 3×3 elements together yields one slice of the new array H , shown on the right of the figure.
For the blue pixel, only 4 elements contribute to the total intensity, so the remaining part of the slice in H remains zero. The backprojection array has the same size and dimensionality as the PSF. Note, however, that the physical meaning of the dimensions has been swapped: The new array now holds 3×3 slices for individual pixels in ( , ) image space, each composed of 3×3 elements attributed to ( , ) object space coordinates. This description formulates a basic recipe for deriving H from H: The forward projection of the individual PSF patterns is actually carried out, storing their contributions to the total intensity of individual image pixels. This approach is realized in the code implementations of [6,8], and requires numerous nested loops with 2D convolutions of rotated PSF slices. Instead, this paper follows a completely different path and proposes an algorithm which exploits the distinct relation of the elements' positions within the arrays H and H . It is based on the fact that the flow of light along rays through the linear optical system can be reversed, which is an expression of Helmholtz' reciprocity principle [20], without involving any reflections in this simple case. It becomes obvious from of Fig. 6 that the transition between the two arrays does not involve any summations, so that H can be found by purely rearranging the elements of H. A step into this direction has been taken by L et al. [9], albeit in a different context, which requires some additional background information. This publication seeks to transfer the PSF (and the recorded light field) into phase-space, i.e. a space/frequency representation. In their light field microscope, the image sensor is placed at the Fourier plane of the MLA, and the pixels within a microimage sample the 2D spatial frequencies of the observed scene. With regard to Fig. 3, this means that here the outer coordinate system defines ( , )-spatial positions, while the inner system is in ( , ) frequency coordinates. In order to transfer H into phase-space, this representation has to be turned inward-out, ordering the spatial information according to its frequencies. This reversal of the physical meaning of the array dimensions is also the result of the transition from H to H , see Fig. 6. It is therefore not surprising to realize that the results of the H computation of e.g. [6,8] are identical to the phase-space transformation used in [9]. The respective code is here based on a restructuring of the original array. It requires, however, an intermediate step where it is expanded into a 6-dimensional form, with a loss in computational efficiency. Before going into the details of the proposed new approach, the following section clarifies why a quick computation of H , though not directly affecting the result of the volume reconstruction, is highly beneficial for light field deconvolution.

Relevance of quick computation of the backprojection array H
It might be argued that the computationally expensive step in light field deconvolution is the iterative volume reconstruction, and deriving the backprojection array is insignificant in comparison. However, modern plenoptic cameras feature a high pixel count, MLAs with hexagonal arrangement and sometimes several lenslet types, which enlarges the necessary elementary cells, see Fig. 4. Using data from such instruments with the published codes leads to very long computing times; in practice, when working with the code from [6] and an R29 multifocus camera from Raytrix, transformation of the light field PSF took 7 h, the same time required for a complete volume reconstruction with 8 iterations. As H has to be computed only once and then can be used for all successive deconvolutions, this may still be regarded as a minor flaw. But there is a number of reasons why a quick computation of the backprojection array is highly beneficial. Especially in cases where the light field PSF is determined experimentally, using a microscope or a photographic setup, its quality has to be assessed after acquisition. Here analysis of H yields valuable hints: Fig. 7 shows a slice of it for an incorrect setting of the physical boundaries and step sizes in object space. The slice reveals clear artifacts and the need for readjustment, which would not have been directly obvious from the patterns stored in H. Further quality enhancement of light field deconvolution could be gained by incorporating physical models into the reconstruction process, that account, e.g., for light refraction due to density gradients within the volume. This would also affect the point spread functions and hence require an update of H -and consequently of H -during each iteration. The benefit of an efficient algorithm is here obvious. Along the same line, we speculate that further developments could lead to blind deconvolution techniques for light field data, where PSF and volume are estimated simultaneously from recorded measurements. Such algorithms are often variants of traditional non-blind techniques, e.g. of the Richardson-Lucy scheme that is common to light field deconvolution approaches [21]. Again, iteratively updating the PSF would also mean computing new backprojection arrays, which requires rapid solutions. In the following, an efficient algorithm for computing the backprojection array H from the light field PSF H is proposed, which is also capable of handling non-symmetric arrays, important for the case of MLAs with hexagonal arrangement. Using sample data, the computational performance of the new algorithm is then benchmarked against previously published codes.

Algorithm
Recalling the notation of Fig. 1, the image plane ( , ) is discretized with × = pixels and the object space ( , , ) is defined by × × = voxels, with the z-axis aligned to the optical axis. In this example, two single voxels at positions ( , , ) are projected onto the image plane, where they are recorded as two-dimensional ( , )-patterns and form respective slices of the PSF array H( , , , , ).
To be more explicit about the details of the image formation process, Fig. 8 again outlines the interplay of the different spaces, for the sake of simplicity only in one spatial dimension ( ), one pixel coordinate ( ) and a single -plane. On the left of the figure, the PSF H has 3 pixels ( = 3) and 3 spatial positions A, B and C ( = 3) colored in red, blue and green. Voxels in object space g (circles) are projected (full lines) according to their respective slices of H to yield the image f. Here dashed boxes represent sensor pixels, where the contributions (filled rectangles) of the different voxels are summed up. For each pixel, H defines the influence of the various voxels, or alternatively, a pixel backprojection into object space. This is sketched by dashed arrows for the middle pixel. In most cases, the image of a single point in object space will have a higher number of pixels in the -and -dimension than there are slices of H in the -and -direction. This means > and/or > and has important implications, which is illustrated in Fig. 8 on the right. Here the PSF H has = 5 and = 3, and this array defines the contributions of 3 points in object space (full circles) being projected (continuous lines) onto each 5 pixels on the camera sensor (full squares), where they sum up to form the image (dashed boxes). To reverse this process, we need to model 3 pixels being backprojected onto 5 object points, and store this in H . This means to include two neighboring points (open circles) and consider their contribution to the image. As the PSF patterns are periodically repeating due to the regularly arranged microlenses, we add shifted slices of H (open squares). Here elements of the lowest slice A, colored in red, reappear on the upper end of H , green elements of slice C at the lower end. With two spatial dimensions and (and, correspondingly, two image dimensions and ), this effect can be modeled by adding slices around H in such a way that both in rows and columns the different slices are repeating in a circular fashion. This is sketched in Fig. 9 for a sample matrix H having = = 3. The core slices in bold letters within the inner rectangle are completed by additional contributions on the sides. In total, ( − ) and ( − ) slices are added in The found relations are at the core of the new algorithm and are going to be detailed in the following. The -planes of the arrays are independent, so no cross-talk has to be considered and they can be computed separately.
With the loop variables and and the auxiliary variables the first step is to derive a spatial position under consideration: Here and denote the and operation, respectively, and the procedure is looped for = 1.. + ( − ), which corresponds to adding ( − ) extra slices. This has to be carried out likewise for the other spatial dimension, resulting in a value using , and . Based on and from Eq. (4), corresponding positions within the array H can be computed as within the limits 0 < ≤ and 0 < ≤ . The relation between dimensions / and / can be found according to with a Matlab version provided in [22]. Separate treatment of and in Eq. (5) allows for non-symmetric cases with ≠ , e.g. with hexagonal MLAs, see Fig. 4. Both new and previous algorithms require the spatial dimensions and to be odd-numbered. A result of the transformation is presented in Fig. 11. It shows two slices of the array H , computed from the experimental data given exemplarily in Fig. 3. These slices define backprojections of single pixels into object space, at identical depth planes, but different lateral ( , ) pixel positions. The off-center slice on the right of the figure clearly shows line patterns due to diffraction at the main lens aperture. This transformation is reversible, so that processing H yields the original array H. As a side note, if all slices of H within one depth plane are summed up, the result is identical to the image given in Fig. 5, except for a rotation by 180 • (or flipping both left/right and up/down), which is due to reversing the projection direction.

Performance
The performance of the new algorithm is tested with sample arrays H of various sizes and is benchmarked against the published codes. In [8], computation of the backprojection array is done with a code that takes advantage from sparse matrices and symmetries within the PSFs. All other algorithms, including the one proposed in this work, do not distinguish between different input data types. Sample PSFs H were therefore generated, based on waveoptics and a microscope system, with routines from [8]. This ensures favorable conditions for their code without compromising the others. In all test cases, PSFs are based on an MLA with a rectangular grid and a single lens type. All algorithms are written in M (version 2019a). Computed backprojection arrays H of the new algorithm and of the codes in [6,8] are completely identical for all samples. Results from [9], designed for phase-space transformation, differ slightly due to cropping the slices by some pixels at the borders, but nevertheless serve as a comparison.
Backprojection arrays H were computed on a desktop PC having an Intel i7-6700K CPU at 4 GHz and 32 GB of memory. Required computation times for all codes, measured in seconds, are given in columns 4, 5, 6 and 7 of Table 1. For all tested sizes, the new algorithm is considerably faster. Especially for high values of and , the algorithms of [6,8] are very slow and the presented new procedure allows a significant acceleration. Column 8 lists the achieved speedup, relative to the best-performing code of [9]. Modern commercial plenoptic cameras feature high resolution image sensors with a high number of pixels under each microlens. As an example, an R29 by Raytrix has 31x31 pixel micro images in an hexagonal arrangement and features three different types of microlenses, a layout sketched on the right of Fig. 4. The representative region, the elementary cell, for such a lens pattern is indicated as a dashed rectangle, and here requires to consider 95x55 positions in the and dimension. The case in the last row of Table 1 is an example for a PSF array H acquired experimentally by a photographic R29 camera, with two slices shown in Fig. 3 and parts of H in Fig.11. Here and are not equal, so that this case cannot be treated by the other algorithms. While the code from [8] in general is capable of handling hexagonal, multi-lens data, PSFs from the different lens types have to be processed separately, which is complicated using experimental data. The new algorithm computes this real-world test case in less than 6 seconds.

Conclusion
This paper has discussed the significance and structure of the array H , which defines a backprojection of an image pixel into object space and forms a key requirement for light field deconvolution methods. This motivated the development of an efficient method to derive H from the shift-variant PSF H, which is commonly given as a 5-dimensional array. It was shown that the positions of individual elements within the two arrays are tied by unique relations, which can be exploited to efficiently compute H from H and vice versa. A new algorithm, based on these findings, was presented, with favorably short computation times compared to other procedures that have been published as part of deconvolution codes. It handles arbitrary PSFs, independent of the lenslet arrangement within the camera's microlens array. A quick calculation of H is beneficial for assessing experimentally acquired PSFs, which often require several adjustments. The general trend towards higher pixel resolutions of digital imaging sensors also holds for plenoptic light field cameras, with commercial devices available in the range of over 100 megapixels. Associated very large PSF matrices of such future systems can be efficiently processed with the algorithm derived in the present work. This could enhance the iterative volume reconstruction by implementing physical models, which account for effects such as refraction due to density gradients, and require to update PSF and backprojection array frequently.
Funding. This work was funded by the German Research Foundation (DFG) under grant No. Lo1772/4-1.