On the Stability of Reconstruction of Irregularly Sampled Diffraction Fields

This paper addresses the problem of reconstruction of a monochromatic light field from data points, irregularly distributed within a volume of interest. Such setting is relevant for a wide range of three-dimensional display and beam shaping applications, which deal with physically inconsistent data. Two finite-dimensional models of monochromatic light fields are used to state the reconstruction problem as regularized matrix inversion. The Tikhonov method, implemented by the iterative algorithm of conjugate gradients, is used for regularization. Estimates of the model dimensionality are related to the number of degrees of freedom of the light field as to show how to control the data redundancy. Experiments demonstrate that various data point distributions lead to ill-poseness and that regularized inversion is able to compensate for the data point inconsistencies with good numerical performance.


Introduction
Many optical applications require a generation and control of light fields.Digital processing by computers is an attractive way of implementing such operations as it overcomes possible physical limitations of analog devices.However, signal processing has to be suitably coupled with the otherwise naturally continuous optical signals.This coupling requires effective discrete representation of the continuous functions, associated with the light fields.From one point of view, the discrete representation of a light field should preserve the degrees of freedom of the continuous model which describes the physical properties of the field.From another point of view, the discrete representation should admit the requirements on the light field which are imposed by the application.Three-dimensional (3D) imaging deals with the reconstruction of captured optical signals by digital means as CCD arrays and recreation of synthetic or computer generated, 3D data by holographic means [1].Light beam shaping requires the reconstruction or synthesis of light beams which maintain certain properties along their propagation [2].
This paper addresses the problem of monochromatic light field reconstruction from ensemble of data samples with free positions, distributed within a volumetric region of interest.Setting the light field specification at a nonuniform irregular grid of points is rather general and can be utilized in a wide range of applications.For example, computer generated holography [3] reconstructs a light field to approximate a synthetic or a captured object, given its 3D abstract model-a point cloud, mesh, NURBS, and so forth [4].The data provided by such 3D models can be hosted conveniently by an irregular grid.Limited aperture and resolution in digital holography are remedied by multiple CCD recordings in a single or different planes [5,6].These are prone to distortions originating from not precisely known CCD measurement locations which can be handled by an irregular grid of points.Reconstruction of a light field in terms of light field generators is required to drive properly a light modulation system for an eventual light field synthesis.For example, a color holographic 3D display system drives three spatial light modulators-one for each red, green and blue color channel [7].
Another advantage of specifying the application constraints as an ensemble of irregularly distributed samples is the control on the input data redundancy.According to the Rayleigh-Sommerfeld diffraction integral, the light field on a certain plane is sufficient to compute the wave field at any point within the volume of interest [8].Therefore, a regular grid of samples at a certain transversal plane should be sufficient to reconstruct the whole field.Gori [9] and Onural [10] have derived conditions for the sampling rate when the support of the field at any Fresnel plane is limited, considering the cases when the diffraction field originates from strictly band limited or strictly space limited function.In the general case of a function which is essentially limited both in space and frequency domain, the amount of needed samples should not be higher than the degrees of freedom N of the field, measured by the finite spacebandwidth product.In this case the conditions from [9,10] require redundant amount of samples, as the spatial spread of the field extends upon propagation along the longitudinal direction z, while the frequency support remains the same.An adhoc approach as in [11] and a strict derivation as in [12,13] use information about local bandwidth to reach the number of degrees of freedom by the amount of needed samples at any transversal plane along z.Early works [3,14] use some of these results to specify the object on a plane, loosing the volumetric opportunity.Some recent works commonly specify the volume of interest on an ensemble of equally spaced planes, at a regular grid of points at each plane, often oversampled [2,[15][16][17][18][19].However, such a 3D grid is still redundant with respect to the number of degrees of freedom of the diffraction field.Moreover, a fixed 3D grid might not fit well to an application-driven data point distribution.In this sense, an irregular grid of points not only fits any potential application but also opens up opportunities for volumetric specifications and data redundancy control.
In any application, the specifications for the light field are not guaranteed to be physically consistent with a light wave.Digital recording of holograms include CCD noise, misalignment, finite aperture and resolution.Light beam shaping and computer generated holography specify complex structures with properties, unlikely to be exhibited by a light wave.The mainstream of articles on light field synthesis treats the data specification inconsistencies by various iterative projection methods, derived on the base of the Gerchberg-Saxton algorithm [20].Volumetric specifications usually define one constraint set per plane, assign weights and tolerances to the sets and use parallel or serial projections to reach an optimized solution [2,6,18,19,21,22].Other methods define the constraint sets on Ewald's surface [17] or fractional Fourier domain [23].Alternative approaches to treat the data inconsistency are based on genetic algorithms with application to beam shaping [24,25].In our work we employ two finite-dimensional models for diffraction which represent a light field by a linear combination of generating functions [26,27].In these models, the inconsistencies of the specified data points can be represented by an additive Gaussian term.While being simple, such a representation is quite general and encompasses a wide range of distortions.
The light field reconstruction problem can be cast as an inverse problem [28] to find the coefficients of the finitedimensional models from the specified data points [26,27,29].The focus of this paper is on the characteristics of the input data specification which influence the stability of the finite-dimensional models for light field reconstruction.The spectrum of each model is analyzed in various scenarios depending on the size of the volume of interest, clustering of the data points and amount of input data, as these are important for practical applications.As far as the mathematical approach is concerned, the Tikhonov regularization [28] is used with both models to provide an approximate solution which fits the specified data points.The singular value decomposition-(SVD-) based implementation of Tikhonov regularization is used to illustrate and analyze the limits of the reconstruction.The regularization parameters are determined by Morozov's discrepancy principle.The matrix-free Tikhonov regularization is suggested as practical approach for reconstruction.It iteratively finds a regularized solution based on the method of conjugate gradients (CGs) [30].The convergence rate and final reconstruction error for various scenarios are presented to assess this approach.
The paper is organized as follows.Section 2 introduces the basics of diffraction and notations related to the considered light fields.Section 3 presents the finite-dimensional models used to represent the field, and discusses their properties with focus on dimensionality issues.Section 4 states the observation model in a matrix form, based on the irregularly distributed inconsistent data samples.Regularized solutions used for approximating the input data are presented in Section 5. Section 6 illustrates the performance of the suggested reconstruction methods in a number of guided experiments.

Basics of Diffraction
Consider a light field u(x, y, z), generated by a monochromatic light wave, which propagates in a linear, isotropic and homogeneous media.Under such conditions, the spatial distribution of the complex amplitude u satisfies the homogeneous Helmholtz wave equation [8]: where k = 2π/λ is the wave number of the monochromatic light.Such waves, emerging from an optical system, satisfy the Sommerfeld radiation condition and are accurately described by the Rayleigh-Sommerfeld diffraction integral [8].It relates the light field u(x, y, z) at any point (x, y, z) to that on a "reference" plane at z = 0 in a linear and shift-invariant relationship.However, the practical use of this relationship for computations is limited as it requires very high sampling rate [31].Therefore, frequency domain alternatives of the Rayleigh-Sommerfeld diffraction integral are preferred [2,22,31].
For the sake of simplicity, the discussion is restricted to one transverse dimension (x) only, resulting in a twodimensional (2D) scalar function u(x, z) describing the light field.Note that the interesting dimension for the scope the paper is the longitudinal one (z) as it keeps the volumetric properties of the light field.Generalizations to the three dimensional case are straightforward and mentioned whenever necessary.The function u(x, z) can be represented by its 2D Fourier transform A(k x , k z ): where k x and k z are the spatial frequencies along x and z directions, respectively.Applying the Fourier transform on (1) and using the derivative property of the Fourier transform, one ends up with the relation ) can be nonzero only on the circle with radius k = 2π/λ: (3)

Finite-Dimensional Models of Diffraction
3.1.Bessel-Fourier (B-F) Generators.As the spatial frequencies k x and k z are limited only to a circle, the Fourier transform A(k x , k z ) is in fact one-dimensional function: where the frequencies are expressed in polar form as k x = k sin θ, k z = k cos θ with k = 2π/λ and θ ∈ [0, 2π).Upon a change of the spatial coordinates x and z in polar form as x = r sin φ and z = r cos φ, the Fourier transform in (2) becomes u r, φ = 2π λ 2π 0 C(θ)e j(2π/λ)r cos(θ−φ) dθ. ( The function C(θ) is 2π-periodic.Consequently, it can be described by the complex Fourier series as C(θ) = m c m e jmθ .Inserting this into (5) and changing the order of summation and integration, the field can be expressed as where J m (t) is the mth order Bessel function of the first kind [32] which solves the remaining integral under the summation.Equation ( 6) assumes that the Fourier series representation of C(θ) has been truncated to M nonzero coefficients c m / = 0, m = − M/2 , . . ., (M − 1)/2 to arrive at finite representation.Detailed derivation of the model can be found in [26].This model describes the field continuously at any spatial point as a superposition of a discrete set of basis functions ψ m (r, φ) = e jm(φ+(π/2)) J m (kr), which are mutually orthogonal and separable.Hence, the discrete set of scaling coefficients c m describes completely the continuous field.Yet, no explicit discretization was done during the derivation.
The model of ( 6) can be generalized to three dimensions [33].In 3D, the Fourier transform A(k x , k y , k z ) is supported on a sphere with the same radius k = 2π/λ and therefore it is two-dimensional.Such a function can be expressed as a superposition of spherical harmonics, instead of the cylindrical harmonics e jm(φ+π/2) .In radial direction the Bessel functions J m (kr) are substituted by the spherical Bessel functions of the first kind.

Fourier Generators.
Another frequency domain representation of the diffraction field u(x, z) is based on a decomposition in plane waves e j(z x z e jkxx dk x , (7) where a(k x ) is the 1D Fourier transform of the field u(x, 0) restricted only to the initial line at z = 0.The derivation of this integral is based on solving the Helmholtz equation for the 1D Fourier transform a z (k x ) of the field u(x, z) restricted to a line at a distance z, given a(k x ) as initial condition.
The plane waves e j(z x +xkx) propagate only in positive direction z ≥ 0 as the frequencies k z = k 2 − k 2 x assume only positive sign.This corresponds to taking only half of the circle k 2 z + k 2 x = k 2 for k z ≥ 0. As k x is limited within [−2π/λ, 2π/λ], the function u(x, 0) is band limited.It can be further assumed to be essentially space limited or, more precisely, to have a finite space-bandwidth product.Such a function can be periodized with period equal to the transversal extent T of a spatial region of interest.This is equivalent to discretization of a(k x ).A periodic and band limited function is equivalently represented by a finite number (M) of discrete spectral components at values of the frequency k x = 2πm/T, m = − M/2 , . . ., (M−1)/2 .Using these values of k x in ( 7) leads to a finite-dimensional model of diffraction: where a m = a(2πm/T) are the coefficients of the Fourier series expansion of u(x, 0).This discrete, Fourier generatorsbased, model describes the field continuously at any spatial point as a superposition of the basis functions ϕ m (x, z) = e j(2π/T)z √ (T 2 /λ 2 )−m 2 e j(2π/T)mx .The generalization of the model in ( 8) is straightforward.Another transversal dimension y brings an extra term in the plane waves e j(z x +xkx+yky ) [8].Now the 2D function u(x, y, 0) must be assumed space limited also along y, as discretization of k y corresponds to periodization along y.

Dimensionality of the Models and Field Characteristics.
The derived models involve a finite number of M basis functions in a linear combination to build up a light field.The amount of the nonzero weighting coefficients determine the dimensionality of any forward or inverse problem related with computation or reconstruction of a light field.Therefore, it is important to relate this amount to the physical properties of the approximating field.The specified region of interest and target detail level can be naturally related to the spatial and frequency content of the field, that is, to the number of degrees of freedom of the field.As the field can be completely described by the function on the initial line u(x, 0), its number of degrees of freedom is equal to the degrees of freedom of u(x, 0).The number of degrees of freedom N is a set of N numbers which describe it completely.In terms of the Wigner distribution, this is the area in the time-frequency plane under which the Wigner distribution of the function is essentially nonzero [34].This can be measured by the space-bandwidth product between the spatial Δx and frequency Δk x extent: where the factor 2π renormalizes the angular frequency k x to an ordinary one.The number N can be related to the number of independent Gabor atoms with unit space-bandwidth product in the Gabor representation of u(x, 0).Alternatively, if u(x, 0) is sampled uniformly along x according to the Nyquist rate 2π/Δk x , then N becomes the number of essentially nonzero samples from which the continuous function u(x, 0) can be reconstructed.In practice, N can often be estimated from the application at hand.For example, it can be related to the physical properties of a sensing device.For synthetic data and other applications, reasonable assumptions about the desired detail level of the target can be used.

Fourier Generators.
The amount of required coefficients M F for this model is directly related with the bandwidth of u(x, 0): Equation (10) suggests that M F is directly proportional to N. This result is expected as the Fourier generators discretize the frequency band of u(x, 0).The proportionality coefficient T/Δx defines the excess which needs to be taken.The period T must be larger than Δx so that the periodic replicas of u(x, 0) do not overlap.However, choosing T close to Δx does not guarantee that an overlap will not occur on a line at further distance z.Discrete frequency axis k x defines discrete Fourier spectrum of u(x, z) at any line z, keeping u(x, z) periodic with the same period T. On the other hand, spatially limited pattern u(x, 0) tends to spread its transversal support when propagated along z.Therefore, given maximal distance z max according to the specified region of interest, one must ensure that T > Δ zmax x, where Δ zmax x is the support of u(x, z max ).The relative increase of this spatial support Δ zmax x/Δx is smaller (closer to 1) for larger Δx at fixed distance z max .

Bessel-Fourier Generators.
The coefficients c m in the Bessel-Fourier model are the Fourier series coefficients of C(θ).Their amount M B can be determined from the frequency support of C(θ).Direct comparison of ( 2) and ( 7) yields Hence, the frequency support of a(k sin θ) coincides with the frequency support of C(θ) and can be used to estimate the number M B .The duality property of the Fourier transform FT{a(k x )} = u(−2πx, 0) can be used to obtain the frequency support of a(k x ) as Δx/2π.However, the frequency support of a(k sin θ) with respect to θ is not the same, even though it is related to Δx/2π.The highest frequency in the Fourier transform of a(k x ) is Δx/4π-the same as the frequency of the harmonic cos((Δx/2)k x ).Therefore, the frequency support of cos((Δx/2)k sin θ) = cos((πΔx/λ) sin θ) can be used to estimate the frequency support of a(k sin θ).In communication theory, such a harmonic function is recognized as a special case of frequency modulated signal cos(2π Its frequency support is approximated as 2( f m + Δ f ) according to the Carson's rule [35].The frequency support of cos((πΔx/λ) sin θ) and a(k sin θ) and C(θ),respectively, is estimated as 2((1/2π) + (Δx/2λ)) ≈ Δx/λ.Hence, the dimensionality of the Bessel-Fourier model required to describe a field with transversal support Δx on the reference line is: This result suggests that the Bessel-Fourier model requires small amount of generators when the spatial support of the field is comparable to the wavelength.From another point of view, if N is assumed to be fixed, then the excess in M B compared to N is determined by the ratio between 2π/Δk x and λ.In fact, 2π/Δk x is the size of the finest detail structure in u(x, 0).

Irregularly Sampled Diffraction Fields
An irregular grid is specified as the set of s sampling points {(x i , z i )} s i=1 , or, in polar coordinates {(r i , φ i )} s i=1 , with the correspondence x i = r i sin φ i and z i = r i cos φ i .The sampled field is specified as the sample values u(x i , z i ), i = 1, . . ., s.The light field reconstruction problem is to find the unknown field-generating coefficients c m or a m , given the irregularly distributed samples u(x i , z i ), i = 1, . . ., s. Equations ( 6) and ( 8) can be written for each point in the irregular sampling set to obtain a m e j(2π/T) for i = 1, . . ., s.Each of these sets of equations forms a linear system for the M unknown coefficients c m or a m and can be expressed in a matrix form: where ] T is the unknown vector of the field generating coefficients and the vector of given samples is u = [u(x 1 , z 1 ), u(x 2 , z 2 ), . . ., u(x s , z s )] T .A is the reconstruction matrix which has two different forms depending on the discrete model.For the Bessel-Fourier model A is and for the Fourier model A is The straightforward approach to solve for the coefficient vector h is to invert the matrix A. Note that the structure of this matrix is determined only by the positions of the known samples.Some field specifications might use more samples to describe the fine details of a scene and less for the uniform regions.Such clustering of the samples causes large condition number of the matrix A and inversion becomes numerically unstable.In addition, there might be noise and/or the scene samples can be inconsistent with the physical models of ( 6) and (8).As the inconsistencies might be of various origin and difficult to predict in the common case, it is reasonable to assume that their nature is also random as any eventual noise appeared upon the scene capture.These two factors can be accounted in a common random term ε: where now y = [y 1 , . . ., y s ] T is the vector of given data samples, h is the unknown coefficient vector, A = J or A = R depending on the model chosen for reconstruction, and ε = [ε 1 , . . ., ε s ] T is a random vector drawn from zero-mean normal distribution ε i ∼ N (0, σ 2 ).Upon direct inversion, the high condition number of A causes strong amplification of the random term and it dominates over the reconstructed coefficients h.

Regularized Reconstruction
A linear system can be solved by taking the pseudoinverse with the help of SVD [36].SVD decomposes an s × M matrix as A = UDV T with U and V containing the left and right singular vectors of A, respectively.D = diag{d 1 , . . ., d n }, n = min(s, M), is an s×M diagonal matrix with the value ordered singular values of A. The Moore-Penrose pseudoinverse A + can be used to find a solution h, provided the desired data vector y as h = A + y = VD + U T y, where If the SVD pseudoinverse is applied to the left side of the noisy observation vector y = Ah + ε, the noise norm in the resultant reconstruction will be boosted by a factor of d −1  n .In the case when the condition number of A is large, d −1 n is large as well, and the noise will be dominant in the reconstruction.A common remedy is to apply Tikhonov regularization, that is, to introduce a penalty term δ into the inversion matrix D + δ = diag{d 1 /(d 2 1 +δ), . . ., d n /(d 2 n +δ)} [28].This is equivalent to solving the minimization problem that is, to minimizing the norms of both residual and solution.Tikhonov regularization balances between the contradictory requirements for small residual and small norm of the solution by the penalty term δ [28].The minimum norm requirement ensures smoothness of the solution an thus robustness to noise, while minimizing the residual ensures that the solution is close to the target.Optimal value of δ can be determined automatically by the Morozov's discrepancy principle [28].
The SVD form of the Tikhonov regularized solution is very convenient for analysis purposes.However, finding the SVD of a matrix is of cubic complexity [36].A more computationally attractive approach is based on iterative matrix solvers [28].The minimizer h of the Tikhonov functional in ( 18) is equivalent to the solution of the linear system [28]: The matrix H = A T A + δI is symmetric and positive definite.Therefore, it can be inverted with an iterative method, which builds the solution step by step, updating the solution vector each time until a desired accuracy is achieved [30].In many cases the iterations converge fast enough to ensure lower complexity than the SVD-based Tikhonov solution.
The conjugate gradient method is one of the most rapidly convergent and numerically stable algorithms [30].It iterates as follows: (1) initialize b = A T y, h arbitrary, r 0 = b − H h and d 0 = r 0 . ( CG method updates the solution at each iteration with a small portion β i along the search direction d i (step (2c)).
The search directions are conjugate to each other, which guarantees convergence in at most M iterations.Practically, the algorithm converges to a predefined accuracy in much less iterations i max M. The most computationally expensive operation at each iteration is the matrix-vector product at step (2a) which requires M 2 operations.The total complexity of the algorithm is O(i max M 2 )-lower than O(M 3 ) required by the SVD-based Tikhonov solution.

Experiments and Results
The light field reconstruction approaches presented in this paper are evaluated by controlled experiments on simulated data.The experiments are described in two subsections.The first subsection describes the details of input data simulation and the criterion for reconstruction performance evaluation, common to all experiments.The second subsection presents experimental results and discusses the benefits of the regularized reconstruction in various scenarios.

Data Simulation and Performance Evaluation.
All experiments are carried out on synthetic data.The data is simulated using the observation model of (17), where the matrix A is computed according to (15) for the B-F model or according to (16) for the Fourier model.The matrices relate M generating coefficients with s data samples, pseudorandomly distributed within a volume of interest.
The dimensionality M of the reconstruction matrix must be the same for both models so that the reconstruction results are comparable.However, this dimensionality depends on different factors for the different models (cf.Section 3.3).For both models, M = 256 nonzero coefficients are able to produce a Gaussian beam inside a square region of interest with size T × z max = 0.1 mm × 0.1 mm.The size of the region of interest is selected so that the adjacent replicas of the field produced by the Fourier generators do not overlap at distance z max comparable to T. Figure 1 depicts the magnitude of the fields produced by the B-F generators and the Fourier generators by showing the values on a uniform 256 × 256 grid which spans the whole region of interest.
The samples of the field used for its reconstruction are simulated by (17) on s randomly selected positions within the volume of interest.The amount of points s is chosen to be greater than the amount of nonzero coefficients M. The underdetermined case s < M is not considered interesting as the minimal norm solution might diverge substantially from the true underlying coefficient vector [27,37].In some experiments, part of the specified data points are selected to form clusters with size much smaller than the region of interest.The samples are randomly distributed according to uniform distribution both inside and outside of the clusters.The size and position of the clusters and the amount of points inside and outside the clusters define different sample density in different subdomains of the region of interest.Varying these parameters helps to investigate the effect of variable sample density within the light field support.The influence of these parameters on the ill-posedness of the reconstruction problem is examined during the first group of experiments.Based on the results, representative scenarios are selected to test the reconstruction approaches.In a fair experiment, the reconstruction performance must be compared with the inconsistency of the input data.Hence, the inconsistency level σ and the reconstruction error e must be measured according to the same quantity.Often one model is used to simulate the samples and another to reconstruct the light field, so that an inverse crime does not take place when experimenting on simulated data [28].The models are characterized by sets of generating coefficients which cannot be compared directly.Instead, the sets of coefficients are used to compute the reconstructed light field within the whole region of interest on a uniform, rectangular grid Γ = {lξ, nζ} M l,n=0 , with ξ = T/M and ζ = z max /M being the sampling steps.The error is measured on this uniform grid: Here, F 0 = u 0 (Γ) and F = u(Γ) are vectors formed by the field values computed on the grid Γ out of the simulation and reconstruction models, respectively.An inconsistency level σ, comparable to this error, can be selected as a percentage of the energy of the field F 0 .

Experiments.
The light field reconstruction problem becomes ill-posed and requires regularization when the spectrum of the reconstruction matrix (R or J) is widely spread and contains many small values.The structure of the matrices depends only on the specified data point distribution.A first group of experiments investigates the spectra of these matrices for different data point distribution scenarios.A second group of experiments illustrates the need and benefit of regularization even for some cases of welldistributed and physically consistent data points.Finally, the last group of experiments shows the benefits and limitations of the regularization for data distributions which make the reconstruction problem ill-posed.This group includes also an experiment which demonstrates a potential practical application of the considered reconstruction approach.sizes and specified data point densities.Spatial variations of the density are simulated by forming data point clusters of various size, shape and position.The plots in Figure 2 represent the singular values of J and R for varying size of the region of interest.It has square shape and the amount of irregularly distributed points is selected as s = 1.5M = 384 so that the matrices are overdetermined.No clusters are formed in this experiment.The vertical axis in Figure 2 represents the singular value d i , while the abscissa represents its index i = 1, 2, . . ., n after the singular values have been sorted in a descending order, that is,

Matrix Singular
Both J and R show their most compact spectra when the dimensionality M = 256 fits the size of the region of interest.Driving away from this matching size, the spectrum of both J and R spreads more.The matrix R has well-concentrated singular values, while J is still ill-conditioned for the matching size.Hence, in the absence of clusters, R can still be used for reconstruction without regularization, as the next group of experiments will demonstrate.Irregularly distributed samples in a square region of interest do not have wide angular diversity and the Bessel-Fourier generators do not bring enough information to recover the coefficients from the samples.This causes wide-spread singular values of J.It is possible to demonstrate that the spectra of J become more compact if the volume of interest has wider size along x than along z.
Next experiments investigate how the presence of clusters influence the spectra of J and R. The basic experiment varies the size of a single square cluster, positioned in the center of a square region of interest of size T × z max = 0.1 mm × 0.1 mm.The experiments use s = 384 points, half of which are distributed inside the cluster and the other half outside.The singular value plots in Figure 3 show that the matrix R becomes ill-conditioned similarly to J. A smaller cluster spreads the singular values more for both matrices.However, the spectra of J are less sensitive to the changes of the size of the cluster, while the spectra of R vary considerably.This is another illustration of the fact that the Bessel-Fourier generators produce a well-conditioned reconstruction matrix when the angular diversity of the samples is higher.On the other hand, this also shows that the Fourier model depends more on the transversal and longitudinal diversity of the sample positions, as the cluster size controls the sample density distribution.Some additional experiments vary the position of a fixed-size cluster within the region of interest and the amount of clusters.The amount of clusters does not influence significantly the spectra of both matrices, provided that the sample density within the clusters is kept constant.However, large amount of clusters spread the general distribution of samples, which improves the spectra of R. The spectra of J improve significantly when the cluster position goes closer to the origin, thus providing samples with higher angular diversity.The spectra of R remain invariant to the position of the cluster within the region of interest.
The last experiments of this group investigates how the amount of given points s influences the singular values of J and R. The simulated region of interest has size T × z max = 0.1 mm × 0.1 mm and contains a single cluster of 100 times smaller size, located in the middle.One of the experiments varies the amount of points outside the cluster from 64 to 1024, while the amount of points within the cluster is kept fixed to 192.The other experiment varies the amount of points inside the cluster from 64 to 1024, while the amount of the samples outside the cluster is kept fixed to 192.The minimum total amount of points s is at least 256 in both experiments, so that the matrices are not underdetermined.The singular value plots of J and R show significant influence when the points outside the cluster are varied (Figure 4), while they remain almost invariant to the amount of points within the cluster.More points outside the cluster carry more information for both models and the singular values of the reconstruction matrices become more compact.The Bessel-Fourier generators benefit from new points outside the cluster as they carry new angles, unlike new points within the cluster.The Fourier generators benefit from more points outside the cluster as this increases the overall point density, Advances in Optical Technologies  while more points within the cluster keep the overall density the same.However, this benefit vanishes after the overall density goes above a certain level, which is theoretically sufficient to reconstruct any field for the chosen volume of interest and dimensionality M.

Regularization for Unclustered Input Data Points.
This group of experiments aims at demonstrating that an ill-conditioned reconstruction matrix needs regularization for reconstruction even for well-distributed and physically consistent input data.The benefits and limitations of the Tikhonov regularization method are illustrated for different levels of inconsistency of the input data.The data was simulated inside a square region of interest of size T × z max = 0.1 mm × 0.1 mm with the B-F model.The model uses M = 256 nonzero coefficients.The amount of simulated points is s = 384 so that the matrices R and J are overdetermined.The points do not form a cluster so that R has compact spectrum and J does not (Figure 2).At first, nonregularized reconstructions were done for different levels of inconsistency with errors shown in Table 1.The reconstruction based on the Fourier model shows practically zero error for physically consistent data (σ = 0%), while the reconstruction with B-F generators shows some very small, yet not zero error.For inconsistent data, the Fourier generators are able to reconstruct a field with an error comparable to the inconsistency level.At the same time the  B-F generators model is unable to reconstruct the field even for very small inconsistency (σ = 1%).This is due to the random term amplification by the reciprocal of the small singular values of J, as discussed in Section 5.In practice, the reconstruction from consistent data is done by applying an iterative technique based on the CG method on the reconstruction matrices R and J.The convergence of CG depends on the condition number d 1 /d n of the matrix [30].Hence, CG converges fast for R whose condition number is small and very slow for J which has high condition number (Figure 5).A regularization strategy is needed when the reconstruction matrix is ill-conditioned.Tikhonov regularization changes the singular values of a matrix by the parameter δ (cf.Section 5).Any change in the spectra of the reconstruction matrix impairs the accuracy of an eventual reconstruction.However, the small singular values are increased and the inconsistency random term ε is not amplified much.An optimal value of δ corresponding to a minimal error is found by Morozov's principle.Such value of δ optimizes the condition number of J T J + δI so that the CG-based Tikhonov regularization converges faster with a minimal final error (Figure 6).The final error is very close to the inconsistency level σ.Thus, the regularization of J helps to acheive very similar convergence and error as reconstructing with the well-conditioned R.

Benefits and Limitations of Regularization for Clustered
Data Points.The last group of experiments aim at illustrating the benefits and the limitations of the regularized field reconstruction when the specified data point distribution contains clusters.In such scenarios, both matrices R and J have widely spread spectra.In this context, it is interesting to investigate how the factors affecting the spectra of R and J influence the regularized reconstruction performance.The factors which have the greatest influence on the spectrum are the size of the cluster and the amount of points outside the  The first experiment of this group investigates the performance of the Tikhoniov regularized reconstruction for different levels of inconsistency σ of the input data.The B-F generators are used to simulate the data points inside a square region of interest of size T × z max = 0.1 mm × 0.1 mm.The total amount of points is s = 384 and half of them form a cluster of 100 times smaller size than the region of interest.The cluster is located in the center of the region of interest.The field was reconstructed from the known samples by the SVD-based Tikhonov regularization on the matrix R for different values of the regularization parameter δ.The plots of the reconstruction error e versus δ are shown in Figure 7.The plots show that the Tikhonov regularization is able to compensate completely for the clusterization of inconsistent input data.If the samples are physically consistent, only an approximate solution can be achieved with sufficiently low error rate.Figure 8 shows the convergence rate of the CG-based Tikhonov regularized reconstruction when the level of inconsistency of the input data is 5%.The value of the regularization parameter δ was computed by Morozov's principle.The B-F model leads to slightly faster convergence rate than the Fourier model, demonstrating the existence of an inverse crime.Comparison with Figure 6 shows that both models have similar performance compared to the case of nonclustered data.In this case, the Tikhonov regularization is able to compensate completely for the ill-posedness caused by data clustering.
Next, it is interesting to test if decreasing the size of the cluster further affects the performance of the regularized reconstruction as it affects the spectra of the matrices R and J (Figure 3).In this experiment, the data inconsistency level is fixed to σ = 5% and the size of the cluster is varied.All other parameters of the scenario are the same as in the previous  experiment.Figure 9 shows the convergence of the CG-based Tikhonov regularized reconstruction involving the Fourier generators.The optimal regularization parameter δ for each curve is obtained by Morozov's discrepancy principle.The optimal δ value is smaller for smaller cluster sizes, and hence the condition number of the reconstruction matrix A T A + δI increases and the convergence of the CG algorithm becomes slower.
Another factor which influences the spectrum of the reconstruction matrices J and R is the amount of given data samples outside the cluster (Figure 4).It is interesting to check the influence of this amount on the reconstruction error and convergence of the CG-based reconstruction.The Advances in Optical Technologies  B-F generators are used to simulate the data points inside a square region of interest of size T × z max = 0.1 mm × 0.1 mm with level of inconsistency σ = 5%.192 points form a cluster in the middle of the region of interest with size 400 times smaller than this region.The amount of points outside the cluster is varied from 64 to 1024 so that the total amount of points is always at least 256 and J and R are not underdetermined.Figure 10 shows the convergence plots of the CG-based reconstruction done with the Fourier model.The lower reconstruction error clearly indicates the benefit of having more data points outside the cluster, showing lower minimal error rate.The optimal value of the regularization parameter δ increases together with the amount of points.This decreases the condition number of the regularized reconstruction matrix A T A + δI and speeds up the convergence rate.
The last experiment investigates a sample distribution scenario which occurs for high-resolution light field reconstruction from multiple low-resolution and limitedaperture CCD recordings.Such a scenario is investigated to demonstrate the practical applicability of the proposed reconstruction approach.In addition, this experiment reconstructs a light field which originates from a realistic object as a more sophisticated distribution.The reconstruction is done again with one transverse dimension only for the sake of comparability with the other experiments.The object at the initial line is 256-sample vertical strip from the "Lena" image, which is used as a standard test material in image processing.The bandwidth of this object contains 256 nonzero frequency components, which coincide with the coefficients a m of the Fourier generators model.The light emerging from such an object is propagated within a region of size T × z max = 0.2 mm × 0.2 mm which makes the effective bandwidth of the image Δk x = 256(2π/2.10−4 ) ≈ 8042.5 rad/mm.The CCD sensors are considered to have twice lower frequency resolution, that is, sampling step of 2(2π/Δk x ) ≈ 0.0016 mm and size 0.1 mm-twice smaller than the transversal extent T. This setting results in 64 samples per CCD.Samples of the light field are captured by 6 different CCD recordings obtained at 6 different CCD positions on a transversal line, such that the oversampling factor is 1.5.The positions of the CCDs are selected such that the total amount of samples is distributed denser around the origin and sparser towards the side of the line, as shown in Figure 11.In this manner, the sampling approximates a density distribution, optimal in the sense of [11].The higher sample density at the center forms a cluster which spreads the singular values of the reconstruction matrix R, as evident from Figure 11.Complex values at the CCDs sample positions are simulated from the coefficients a m of the Fourier generators model.In practice, complex values can be obtained by, for example, temporal phase shifting for each CCD position.The CCD noise is simulated by adding an inconsistency with level σ = 5% to the simulated CCD sample values.The regularized, CG-based reconstruction shows rapid convergence to an error comparable to the inconsistency level, as illustrated in Figure 12.

Conclusion
This paper has addressed the problem of monochromatic light field reconstruction from irregularly distributed samples with physically inconsistent values.The proposed reconstruction method is based on finite-dimensional modeling of the problem and regularized inversion.Our approach encompasses a wide range of applications in the area of 3D display and beam shaping.The dimensionality of both models can be directly related to the number of degrees of freedom of the field with a certain excess.The excess depends on the ratio between the transversal extent and volume of interest for the Fourier model and on the ratio between the finest detail level and wavelength for the B-F Advances in Optical Technologies model.Uniform sample density defines the reconstruction problem as well posed when described by the Fourier model and as ill-posed when described by the B-F model.For such density, satisfactory reconstruction can be done with amount of points close to the number of degrees of freedom of the field.However, spatial variations of the density define the reconstruction problem as ill-posed for both models.Regularized reconstruction is able to compensate for the illposedness when the sample density variations are not very diverse within the region of interest.The reconstruction error and convergence of the iterative solver are very similar to the well-posed case.Large variations of the sample density increase the reconstruction error above the level of inconsistency of the input sample data.In this sense, increasing redundancy in the amount of samples away from the number of degrees of freedom of the field brings clear benefit.

Figure 1 :
Figure 1: A diffraction field representing a Gaussian beam computed by the B-F generators (a) and by the Fourier generators (b) on an uniform 256 × 256 grid covering a region of interest with size 0.1 mm × 0.1 mm.

AdvancesFigure 2 :Figure 3 :
Figure 2: Singular values of J (solid lines) and R (dashed lines) for different sizes of the region of interest in [m] (number next to each line) and s = 384 given data points which do not form a cluster.

Figure 4 :
Figure 4: Singular values of J (solid lines) and R (dashed lines) for different amounts (shown next to each line) of given samples outside a cluster of size 100 times smaller than the region of interest.

Figure 5 :
Figure 5: Convergence of the CG algorithm for reconstruction from nonclustered and physically consistent field samples, simulated by the B-F generators.The CG algorithm was applied to the nonregularized matrices J and R.

Figure 6 :
Figure 6: Convergence of the CG algorithm for reconstruction from nonclustered field samples, generated by the B-F model.The inconsistency level is 5%.The CG algorithm was applied to the regularized J and unregularized R.

Figure 7 :
Figure 7: Tikhonov regularization reconstruction error e versus the regularization parameter δ for different noise levels σ.The samples are simulated with the B-F generators and the field is reconstructed by the Fourier generators.The samples form a square cluster 100 times smaller than the volume of interest.

Figure 8 :
Figure 8: Convergence of the CG-based Tikhonov regularized reconstruction from samples generated by the B-F model which form a square cluster of size 100 less than the region of interest.The inconsistency level of the samples is 5%.

Figure 9 :
Figure 9: Convergence of the CG-based Tikhonov regularized reconstruction for different sizes of the cluster (as percent of the region of interest).The samples are simulated with the B-F model and σ = 5% inconsistency and the reconstruction is done by the Fourier generators.

Figure 10 :
Figure 10: Convergence of the CG-based Tikhonov regularized reconstruction with the Fourier generators for different amounts of data points outside a cluster (number next to each line).The samples are simulated with the B-F generators and σ = 5% inconsistency.

Figure 11 :Figure 12 :
Figure 11: Samples from CCD recordings at 6 different CCD positions (shown in different colors and heights) on a transversal line (top) spread the singular values of R (bottom).

Table 1 :
Nonregularized reconstruction errors for different noise levels.