Studies on the sparsifying operator in compressive digital holography

: In compressive digital holography, we reconstruct sparse object waveﬁelds from undersampled holograms by solving an (cid:96) 1 -minimization problem. Applying wavelet transformations to the object waveﬁelds produces the necessary sparse representations, but prior work clings to transformations with too few vanishing moments. We put several wavelet transformations belonging to di ﬀ erent wavelet families to the test by evaluating their sparsifying properties, the number of hologram samples that are required to reconstruct the sparse waveﬁelds perfectly, and the robustness of the reconstructions to additive noise and sparsity defects. In particular, we recommend the CDF 9 / 7 and 17 / 11 wavelet transformations, as well as their reverse counter-parts, because they yield su ﬃ ciently sparse representations for most accustomed waveﬁelds in combination with robust reconstructions. These and other recommendations are procured from simulations and are validated using biased, noisy holograms.


Introduction
Compressive Sensing (CS) deals with the reconstruction of sparse signals from sets of independent linear measurements that are incomplete according to the Shannon-Nyquist sampling theorem [1].When we sample analog signals at a rate that exceeds the Nyquist bound, an additional compression step is necessary to represent the signals accurately by a small set of decorrelated coefficients [2,3].Minimizing the number of samples or linear measurements we have to collect in the first place is the purpose of CS.Efficient image recording procedures have emerged in this context [4][5][6][7][8][9][10][11][12] at the expense of an increase in computational complexity because the sparse signals are ideally the unique solutions of convex optimization problems that are solved using iterative algorithms.
Recently, CS has been applied to digital holographic imaging with success [13,14] because of two reasons: (1) wavefields are sparse in the spatial domain [15] or a(n) (orthogonal) wavelet domain [16,17], and (2) the propagation of wavefields in free space is incoherent with a basis of wavelets if the interferometer is well configured [18].Since signal sparsity, sufficient incoherence, and randomly selected measurements are excellent conditions for CS, we can assert that applying CS to digital holography is promising.
The integration of CS in digital holography is called Compressive Digital Holography (CDH) and has lead to remarkable achievements including the reconstruction of partially occluded objects [19], the reconstruction of three-dimensional tomographic non-diffuse [20] and diffuse [21] objects, increased axial resolution for in-line holography [22], reduced scanning efforts in incoherent multiple view projection holography [23], digital holography with a single-pixel camera [24], compressive holographic video [25], and pixel super-resolution for lens-free digital holography [26,27].
The reconstructed object wavefields are subject to an error that is proportional to the sparsifying operator's ability to concentrate most of the signal energy in very few transform coefficients.Since the reconstruction is the result of an iterative algorithm, fast transformations like the discrete cosine transformation, the total variation operator, or a wavelet transformation are desired.In particular, wavelet transformations should be characterized by sufficient vanishing moments, V , so that wavelet coefficients are (close to) 0 for wavelets that overlap with parts of the signal that can be approximated by (complex-valued) polynomials with a degree of at most V − 1, see Paragraph 6.1.2 in [28].Prior work in CDH solely reports sparsity in orthogonal wavelet bases with at most 2 vanishing moments [14] or casts the reconstruction as a total variation minimization problem [29,30].Especially the Haar wavelet transformation (also called the Daubechies 2 wavelet transformation) with a single vanishing moment is employed very often, whereas state-ofthe-art image compression standards like JPEG 2000 recommend the Cohen-Daubechies-Fauveau (CDF) 9/7 wavelet transformation for lossy compression.This and other wavelet transformations of the same family are best suited to image compression because the corresponding compressed images are the least susceptible to visually disturbing edge distortions [31].In fact, the CDF wavelet transformations constitute just one class of wavelet transformations.We will also take a look at the Daubechies n wavelet transformations for n > 2, the Battle-Lemarié wavelet transformations and the (dual) B-spline 3 wavelet transformations [32].We include the latter because they are the cornerstone for the construction of the Fresnelets, which are promoted for the compression of digital holograms in [16].
The aforementioned works about CDH show only the amplitude of the reconstructed object wavefields, although the phase information is important too.The depth information or threedimensional geometry of the object under test is namely determined from the phase profile.We focus therefore on both the amplitude and phase of the reconstructed object wavefields.
In this paper, we promote wavelet transformations with sufficient vanishing moments and in particular the CDF 9/7 and 17/11 wavelet transformations, as well as their reverse counterparts, as sparsifying operators for object wavefields in a CDH framework.We substantiate our recommendations with both simulations and experimental data, as we show in Code 1 [33].In particular, our extensive studies show that the sparsifying transformations that have been employed in prior works have insufficient vanishing moments.Theoretically, this is not an issue for Piecewise Constant (PWC) object wavefields: we confirm that the Haar wavelet transformation with one vanishing moment is the best sparsifying operator in a numerical simulation, but not in practice because sharp edges are smoothed in biased, noisy experiment recordings.As such, the Haar wavelet transformation should not be the default sparsifying operator in CDH according to us.Selecting wavelet transformations with more vanishing moments results in less blocking artifacts.On the other hand, we find that the (dual) B-spline wavelet transformations are good sparsifying operators for smooth wavefields, but they lead to reconstructions that are susceptible to artifacts we do not observe for the (reverse) CDF 9/7 and 17/11 wavelet transformations (shorthand notation to refer to the four CDF wavelet transformations).As such, we do not promote the (dual) B-spline wavelet transformations for CDH.Finally, we show that for the propomoted wavelet transformations the reconstructions are more accurate than for the discrete cosine transformation or a total variation prior.

Compressive sensing
CS relies on two conditions: (1) sufficient signal sparsity and (2) incoherence between the sensing and sparsifying operators.Let us consider a d-dimensional object wavefield consisting of a linear sparsifying operator of f 0 , then we can model the sensing procedure as where α 0 = Ψ −1 f 0 is the S-sparse representation of f 0 (i.e., α 0 has S non-zero elements) and T is a matrix that represents an index set for identifying the hologram samples that are actually measured.Hologram samples that are no elements of T are thus masked.Due to the homogeneous, compact and regular structure of accustomed objects, the corresponding object wavefields are sparse in the spatial or a transformed domain.Several wavelet transformations are put to the test in this study to find the most sparsifying one for PWC and Non-PWC (NPWC) object wavefields.We note that biorthogonal systems -like the (reverse) CDF 9/7, (reverse) CDF 17/11, and (dual) B-spline 3 wavelet transformations -consist of two biorthogonal bases.One basis is used for decomposing the object wavefield in wavelet coefficients while the other one is used for reconstructing the object wavefield.The two variants are distinguished from one another by the absence or presence of the adjectives 'reverse' or 'dual' in front of the name of the wavelet transformation.For correct implementations we rely on the Wavelab package [34], except for the (dual) B-spline 3 wavelet transformations, which we implemented ourselves based on [32].Furthermore, for the reader's convenience, we list in Table 1 how many vanishing moments the wavelet transformations have.
If Φ and Ψ are orthogonal transformations, the coherence between them is defined as for i, j = 1, 2, • • • , N d .We grasp from Theorem 1 in [1] that the number of randomly selected hologram samples, M, has to be greater than K µS log(N 2d ) for perfectly reconstructing α 0 or f 0 from g T with overwhelming probability by means of an 1 -minimization algorithm.So M decreases as the coherence µ between Φ and Ψ decreases.Unfortunately, the value of the constant K > 0 is unknown and the calculation of µ in Eq. ( 2) is only valid for orthorgonal wavelets and not for biothogonal wavelets.Hence the value of M is hard to predict using Eq. ( 2).We rely therefore on lots of numerical test cases that lead to a relationship between the sparsity of the object wavefield, S, and the required amount of hologram samples, M, such that the sparse object wavefield can be reconstructed perfectly with overwhelming probability.By expressing the probability of success as a function of S and M we calculate so-called phase transition diagrams [35].
Since measurements are never infinitely precise, we solve for the unknown α in the quadratically constrained basis pursuit problem, where the upper bound on the residual's norm, ε, is an estimate of the Euclidean norm of the additive noise.We say that the reconstruction is robust with respect to additive noise and sparsity defects if the Euclidean norm of α 0 − α is bounded.The Restricted Isometry Property (RIP), which is explained shortly here, offers the means to evaluate the reconstruction error.Let us denote TΦΨ by A T .We assume that the non-zero elements of the sparse vector α 0 are known.The columns of A T that correspond with these non-zero elements should constitute a submatrix that is not rank deficient.Otherwise, there is no hope to find α 0 from measurements that are formed by multiplying A T with α 0 .The RIP builds on this reasoning for the case where the exact locations of the non-zero elements are unknown.Given that α 0 is an arbitrary S-sparse signal and c ∈ R + 0 , we say that A T satisfies the RIP with restricted isometry constant In particular, the reconstruction error is a linear combination of the variance of the additive noise and the Euclidean norm of the sparsity defects if δ S < 0.62, see Theorem 6.12 in [36].The reconstruction becomes thus more robust as δ S decreases, but exact values of δ S are unknown (in practice).In the next section, lots of simulated test cases for which we determine δ S give us insight in the robustness for different wavelet transformations while validation is provided by visual inspection in Section 6.

Compressive digital holography in the near field
We select the angular spectrum method as the sensing model for the propagation of waves in free space [37].The wavefields f 0 and g 0 are defined in planes that are perpendicular to the optical axis.The origins of f 0 and g 0 coincide on top of that with the optical axis and they are a distance z apart.Substitution of the expression for the sensing operator Φ in Eq. (1) yields where F {.} is the two-dimensional unitary Fourier transformation, is the coherent transfer function for waves that propagate through free space, and λ is the wavelength.
The discrete form of Eq. ( 5) and its inverse, which is also referred to as the Back Propagation method, are discussed in [38].Both in-line phase-shifting and off-axis methods can be employed to acquire the complex wavefield g 0 [39].
In CDH, the reconstructed object wavefield f 0 is inferred from an undersampled hologram g T .Just as in [40], the samples are selected uniformly at random.This undersampling scheme is also used in single-shot phase imaging with a coded aperture [41][42][43] so the setup we propose and the results of our study on the sparsifying operator are relevant to existing computational imaging methods.Our goal is to identify which wavelet transformations Ψ are preferred such that a good estimate f = Ψ α of f 0 is obtained by solving the convex optimization problem in Eq. ( 3).We selected the SPGL1 [44] solver for this purpose because it comprises a fast algorithm that has already been put to the test in compressive magnetic resonance imaging [8], which is also a complex-valued CS problem.

Undersampling bounds
The coherence µ is a quantity that allows us to estimate roughly the minimum amount of noisefree, unbiased measurements, M, that are required to reconstruct a given S-sparse signal perfectly.
In particular, Theorem 1 in [1] establishes that M decreases as the coherence µ between the orthogonal transformations Φ and Ψ decreases.The angular spectrum method (Φ) and orthogonal wavelet transformations (Ψ) are all orthogonal transformations so it makes sense to compute the coherence µ between them.Since Φ is dependent on λ, z and the image sensor's pixel pitch ∆x, the coherence is also dependent on these three parameters.The laser and camera are off-the-shelf components.Therefore, we set λ = 633 nm (HeNe-laser, JDSU Uniphase 1135P) and ∆x = 9.1 µm, which is 3 times the pixel pitch of the Ximea MD120MU-SY.We observe in Fig. 1 that the coherence decreases as the propagation distance z increases.As a consequence, M decreases as the distance between the object and the camera increases.In particular, the results for the Daubechies 16 wavelet transformation show that µ decreases as z increases, which is in accordance with what has been reported in [18] for CDH with sparsity in the spatial domain.Furthermore, we demonstrate that the trend of the coherence as a function of z is maintained for the CDF 9/7 and reverse CDF 17/11 wavelet transformations.These biorthogonal wavelet transformations are characterized by wavelet bases that are quasi-orthogonal so the coherence is still a meaningful quantity.Unfortunately, we cannot position the camera very far from the object because beyond a certain distance the lateral spatial extent of the propagating wave would be much larger than the image sensors area, which would result in a very small numerical aperture and thus excessive information loss.Moreover, the angular spectrum method would suffer from aliasing artifacts [38].We set therefore z = N (∆x) 2 /λ, which is the maximum propagation distance for which aliasing is avoided.
The coherence is a very common tool in applied CS domains, but there are two issues: (1) the undersampling bound M ≥ K µS log(N 2d ) is only determined up to an unknown constant K ∈ R + 0 , and (2) Φ and Ψ have to be orthonormal transformations [45].As such, we cannot  evaluate the undersampling bound for generic biorthogonal wavelet transformations by means of the coherence.For accurately quantifying the undersampling bounds we rely on phase transition diagrams.
Let us consider a mapping from the phase-space (δ, ρ) ∈ [0, 1] 2 to a probability π, where the undersampling factor δ := M/N d and the oversampling factor ρ := S/M.For each combination of δ = 0.02k and ρ = 0.02 , where k , l ∈ {1, 2, • • • , 50}, we generate a set of R = 50 randomly selected one-dimensional (d = 1) S-sparse signals.The relationship between S and a point in the phase-space is S = δ ρN d .The sparse signals are drawn from a distribution with the following three properties; S out of N = 512 elements are selected uniformly at random without replacement; the values of the selected elements are drawn from a standard circularly-symmetric complex normal distribution while the remaining N − S elements are assigned zero; and the signals are normalized (unit Eucledian norm).Each experiment targets reconstructing an S-sparse signal by solving the optimization problem in Eq. ( 3) with the SPGL1 solver and = 0 (no noise).Notwithstanding that the SPGL1 solver is known to underestimate the undersampling bound slightly [35], it comprises algorithms that are based on a fast spectral projected gradient pursuit for solving large-scale complex-valued reconstruction problems.We stay with the standard solver settings for complex problems, except for 'bptol' and 'optTol' which are lowered to 1 × 10 −6 .
An experiment is labelled 'successful' if the mean squared error between the original S-sparse signal and its reconstruction is smaller than 1 × 10 −12 .We formulate like this the outcome of every experiment as a Bernouilli trial and we assume that all R outcomes are drawn independently from a Bernouilli distribution with success rate π(δ, ρ).By means of maximum likelihood estimation we obtain π where C(δ, ρ) is the number of successful experiments that are associated with the pair (δ, ρ) in the phase-space.
In the last step, the undersampling bounds are derived from the phase transition diagrams by searching in the phase-space for the isoline ρ * (δ) : δ ∈ [0, 1] → ρ ∈ [0, 1] such that π(δ, ρ * (δ)) = 0.5.We show in Fig. 2 the phase transition diagram for CDH with sparsity in a basis of Haar wavelets and the resulting undersampling bound, which is on the verge of the black and white areas.For the other wavelet transformations, we summarize in Table 1 the quantity which measures the Maximum Absolute Deviation (MAD) of the undersampling bound for any wavelet transformation with respect to the undersampling bound for the Haar wavelet transformation.We observe that the MADs are smaller than 0.05, except for the (dual) B-spline wavelet transformations.Their undersampling bounds are therefore depicted in Fig. 2. They are clearly below the undersampling bound for the Haar wavelet transformation, which means that for any fixed sparsity S, more measurements have to be collected when the sparsifying operator is a (dual) B-spline wavelet transformation in comparison with the other wavelet transformations.The Daubechies, Battle-Lemarié, and (reverse) CDF wavelet transformations are associated with undersampling bounds that are quasi-identical and in the proximity of the undersampling bound for the Haar wavelet transformation.The required amount of measurements for these wavelet transformations is thus mainly determined by the degree to which they sparsify a wavefield (i.e., the specific value of S) rather than the locus of the corresponding undersampling bound in the phase space.
The phase transition diagrams specify to which extent the holograms can be undersampled such that the failure rate is sufficiently small.However, they do not explain why the undersampling bounds for the (dual) B-spline wavelet transformations are significantly below the undersampling bounds for the other wavelet transformations.Moreover, we cannot deduce to which extend the solution of optimization problem (3) is robust to real-life wavefields that are only approximately S-sparse and holograms that are corrupted by additive noise.The restricted isometry property is employed here to quantify the reconstruction error.
We note that our formulation in Eq. ( 4) is different from the one in [36].We bound the squared Euclidean norm around a generic constant c rather than 1 to account for eigenvalues of A H T A T that are not concentrated around 1.This is allowed because scaling the (noisy) hologram by β (g T → βg T , → β ) as a consequence of neutral density filters or varying camera gain, leads to the same scaling factor for the coefficients of the reconstructed object wavefield ( α → β α).Let us rewrite Eq. (4) as where H denotes the conjugate transpose operator.We observe that the eigenvalues of A H T A T have to be concentrated around constant c, which is the average of the maximum and minimum eigenvalue, for δ S to be near zero, where we note that c times δ S is equal to the maximum eigenvalue minus c and c minus the minimum eigenvalue.Considering that A T is dependent on a random realization of matrix T, we compute the eigenvalues for 1000 independent realizations.Furthermore, we repeat this procedure for d = 1, S ∈ {32, 64}, M ∈ {128, 256, 512}, and the wavelet transformations that are listed in Table 1.The quartiles of the resulting distributions of δ S are shown in Fig. 3 using box plots, from which we distinguish two observations.Firstly, δ S reaches 0 for orthogonal wavelet transformations when M = N while the CDF and B-spline wavelet transformations score less good with respectively δ S ≈ 0.2 and most δ S > 0.4.Secondly, we find that δ S increases as M decreases.Less holographic measurements lead thus overall to larger reconstruction errors.The gap between orthogonal and CDF wavelet transformations thereby cancels when M < N while the B-spline wavelet transformations still score significantly worse.The restricted isometry constants (values of δ S ) are obviously farthest from 0 for them.We expect therefore larger reconstruction errors for B-spline wavelet transformations in comparison with the other wavelet transformations when they sparsify an object wavefield equally well (with the same sparsity S).This will become apparent in the next sections.

Sparsity-driven CDH study
Since the undersampling bounds and restricted isometry constants are rather insensitive to the type of wavelet transformation, except for the (dual) B-spline wavelet transformations, the reconstruction error is highly dependent on the sparsity of a wavefield in a particular basis.In this section, we investigate the reconstruction error of object wavefields for several wavelet transformations by means of simulations.

Description of the simulations
We consider (partially) occluded and reflective objects that cause a change in the optical path length of up to 10 times the wavelength.The first object is a US Air Force (USAF) 1951 resolution chart that acts as a mask when it is illuminated in transmission.The resulting object wavefield, which is visualized in Fig. 4(a), is PWC.The second object is a lens-shaped reflective surface that leads to the NPWC object wavefield that is depicted in Figs.4(b) and 4(c).Its amplitude profile is circularly symmetric and the sag (maximum thickness) is 6 µm.
The computer-generated object wavefields are instances of f 0 in Eq. ( 5).They are sampled on an N × N square grid (N = 1024, d = 1) with a pixel pitch ∆x of 9.3 µm.These wavefields are illuminated by a plane wave of wavelength λ = 633 nm and propagate over a distance z = N (∆x) 2 /λ such that f 0 is projected onto g 0 as specified by Eq. ( 5).Subsequently, g 0 is undersampled uniformly at random so that the resulting vector of measurements is an instance of g T in optimization problem (3).To complete the instantiation of Eq. ( 3), we note that the sparsifying operator Ψ is any of the wavelet transformations in Table 1.
Optimization problem (3) is solved by means of the SPGL1 solver with the same settings as in Section 4. The resulting estimated wavelet coefficients are then mapped by the sparsifying wavelet transformation to the estimated/reconstructed object wavefield f .For evaluating the reconstruction error we have selected the Signal to Error Ratio (SER), which is defined as where the energy of a discrete signal is the sum of the squared sample values.A reconstructed object wavefield f for which the SER is greater than 150 dB, is labeled as a perfect reconstruction in the remainder of this paper.

Evaluation of the sparsity
To evaluate the sparsity of the computer-generated object wavefields in the wavelet domain, we consider the approximation that consists of the n(p) most-significant/largest wavelet coefficients, where n(p) is the smallest positive integer such that the SER in the spatial domain is greater than p dB. Dividing n(p) by the total number of samples on the grid (N d = N 2 = 1024 2 ) yields a normalized sparsity metric n(p) = n(p)/N 2 for evaluating the sparsity.The results are shown in Table 2 for p = 20 and p = 30.
The PWC USAF resolution chart is most sparse in the Haar wavelets because they have firstly the shortest support size among all wavelets.When an edge or corner of the USAF resolution chart's bars overlaps with the support of a wavelet, the corresponding wavelet coefficient is excited, thereby compromising the sparsity.The other wavelets have larger support sizes than the Haar wavelets so the singularities (edges and corners) can excite coefficients of wavelets that lie in larger neighbourhoods, thereby harming the sparsity even more.Besides, the Haar wavelets are characterized by a single vanishing moment (V = 1).Polynomials of degree V − 1 = 0 (i.e., constant functions) are thus orthogonal to the Haar wavelets, resulting in wavelet coefficients that are 0. Since the transparent areas of the USAF resolution chart do not disturb the incoming plane wave locally, the object wavefield is constant in these areas so that 1 vanishing moment suffices to restrain locally the wavelet coefficients to 0. The trade-off between sparsity and vanishing moments is supported by the values of n(p) for the Daubechies n wavelet transformations in Table 2; as n increases, the wavelets' support sizes increase too and as a consequence the values of n (20) and n(30) increase monotonously.
The opposite trend is observed for the smooth object wavefield due to the reflective lens-shaped surface.No singularities are present in this wavefield so a small support size for the wavelets is of less importance.However, the curvy profile has locally the shape of a polynomial with a degree that is significantly larger than 0. The Haar wavelets fail consequently to sparsify this smooth wavefield.Increasing the number of vanishing moments is required for the wavelets to be orthogonal to polynomials with increasing degrees, thereby mapping the wavelet coefficients to 0 as is explained in Paragraph 6.1.2 of [28].As a result, all wavelet transformations with 4 or more vanishing moments produce values smaller than 0.01 for n (20) 2 are deployed as sparsifying operator.See Data File 2 for underlying values.

The reconstruction error
The quantity of interest is the reconstruction error.Therefore, we solve optimization problem (3) for the two computer-generated wavefields, the wavelet transformations in Table 2, and undersampling factors ranging from 0.05 to 0.95.The resulting reconstruction errors are shown in Fig. 5, where a first inspection reveals that the SER increases as the undersampling factor increases.Simply put, because M is proportional to δ = M/N d (M/N 2 for the two-dimensional wavefields in Fig. 4) and more independent measurements carry more information about the object wavefield.Thanks to this monotonically increasing behaviour we can prune the SER curves to the smallest undersampling factor for which perfect reconstruction is achieved in order to avoid clutter in Fig. 5.
For the Daubechies n wavelet transformations, we find that the SER increases for all undersampling factors as n(p) in Table 2 decreases.The decreasing trend of the reconstruction error as the number of vanishing moments (V = n/2) increases, is apparent in Figs.5(a) and 5(b), where n appears in the markers.Furthermore, we note that the Daubechies 16 and Battle-Lemarié 3 wavelet transformations are characterized by quasi the same undersampling bound (see Table 1), sparsifying behaviour (see Table 2), and restricted isometry constants (see Fig. 3), but the recon-  struction errors differ significantly.We attribute this behaviour to the fact that the reconstruction error is proportional to the restricted isometry constant δ S , but we have no knowledge of the actual proportionality constant, which is dependent on the specific realization of matrix A T and thus also on the selected wavelet transformation and undersampling factor.We learn from this observation that the selected metrics only predict the reconstruction error coarsely.We present the results for the biorthogonal wavelet transformations in Fig. 5(b) to avoid a web of crossing and overlapping lines in Fig. 5.In the class of CDF wavelet transformations we find that the reverse CDF 9/7 wavelet transformation produces the sparsest wavefield representation for the PWC object while both the CDF 17/11 and reverse CDF 17/11 wavelet transformations take that role for the NPWC object.As a result, they lead to the curves with the smallest reconstruction errors (highest SERs).Within this class, we also discern that the reverse CDF 9/7 and reverse CDF 17/11 wavelet transformations produce the smallest reconstruction errors for respectively the PWC and NPWC object.
Finally, we confirm that the (dual) B-spline 3 wavelet transformations are good sparsifying operators, especially for NPWC objects.However, due to the weaker undersampling bounds and less robustness to sparsity defects, they fall short in a CDH problem and lead to the largest reconstruction errors when the undersampling factor is around 0.1.

Convergence rate
When implementing CDH it is also important to keep the computation time as low as possible.Because wavelet transformations can be implemented with several methods like filter banks and lifting schemes, which are characterized by different performances, we have selected the number of iterations to reach convergence as a measure.See Fig. 6 for the results using the most sparsifying operators of the USAF resolution chart and the lens-shaped relfective surface.We observe that less measurements (δ decreases) lead to more iterations.Small δ are desired nevertheless.Therefore, the most sparsifying operator should be selected because we observe in Fig. 6 that this strategy leads to the least iterations in combination with the smallest reconstruction error.

Visualization of some reconstructions
Let us visualize some of the results.We display in Fig. 7 the USAF resolution chart's reconstruction, starting from a hologram that is undersampled 20 times (δ = 0.05).Reconstruction artifacts are salient in the big square between the elements of groups 2 and 3.For the Haar wavelet transformation we observe in Fig. 7(a) that the intensity profile of the big square is flat.On the other hand, the reconstructions for the reverse CDF 9/7 and 17/11 wavelet transformations are affected by ripples in the intensity profiles in Figs.7(b) and 7(c).In other parts of the reconstruction, the artifacts are less outspoken, though still visible.Especially near edges and corners.We selected as sparsifying operator for the smooth lens-shaped surface's object wavefield a wavelet transformation with more than 3 vanishing moments.In Fig. 8, the amplitude and phase of the complex-valued reconstructions (δ = 0.05) for the CDF and Battle Lemarié 3 wavelet transformations are presented.The difference between the CDF 9/7 and reverse CDF 9/7 wavelet transformations is noticeable in Figs.8(b) and 8(c) and is a consequence of the fact that reverse transformation is a significantly better sparsifying transformation (see Table 2).Although they both have 4 vanishing moments, we find significantly more visual distortions in the reconstruction with the CDF 9/7 wavelet transformation as sparsifying operator.Especially the amplitude is highly distorted by rasterization, whereas the circularly-symmetric phase profile is still recognizable.Four vanishing moments are thus not enough in general.The (reverse) CDF 17/11 wavelet transformations with 6 and 8 vanishing moments lead to smaller reconstruction errors according to the SER curves in Fig. 5, but this is barely noticeable in Fig. 8. Finally, the Battle Lemarié wavelet transformation is associated with the smoothest amplitude profile.Indeed, it is the most sparsifying operator according to Table 2 and leads to the highest SER according to Fig. 5.

Experiments
We validate our recommendations on holograms that we acquired in a laboratory setting.A Gaussian laser beam (JDSU Uniphase 1135P; λ = 633 nm) with a beam diameter of 15 mm illuminated three transmissive objects: a USAF resolution chart (flat), a leaf from which the chlorophyl was removed (≈ 4λ thick), and an angle grid (flat).For the results in Subsections 6.1 and 6.2 we positioned at a distance of respectively 4.48 cm and 7 cm from the object a Ximea MD120MU-SY CCD camera (3.1 µm pixel pitch; 4242 × 2830 pixels).Three phase-shifted interferograms of each object were acquired, the borders were cropped, and blocks of 3 by 3 pixels were combined to increase the pixel pitch to 9.3 µm.The resulting holograms have side lengths that are greater than 512 pixels.As a consequence the near field extends to a propagation distance of approximately z = N λ/(∆x) 2 = 7 cm for N = 512, which justifies the position of the camera.Finally, the triple of holograms was processed by the 3-step method from phase-shifting digital holography [46], followed by numerical back propagation to obtain the noisy object wavefields.

Rapid phase variations are difficult to sparsify with wavelet transformations
We recommend wavelet transformations as sparsifying operators for object wavefields, but this does not mean that any wavefield is sparse in the wavelet domain.Being aware of the conditions that harm the sparsity is invaluable.Diffuse objects for example are not sparse in any wavelet domain and are thus excluded from the set of allowed objects.
Let us consider further a tilted plane wave such that the maximum difference between successive pixels is close to π radians.This will introduce pervasive high-frequency content in the hologram, which is suboptimal for Mallat wavelet decompositions: almost all of the signals energy will be relegated to the highpass subbands of the first wavelet level, causing the transformed signal to be insufficiently sparse.These type of signals are more effectively represented with e.g.DCT transforms or by wavelet packet decompositions.
Where the reconstruction error manifests itself in the object wavefield cannot be predicted by existing theorems in the field of compressive sensing.We conducted therefore an experiment with the object wavefield of an array of cantilevers whose depth profiles follow a superlinear inclination.As a consequence, the variation of the phase increases from the bulk to the tips, which is illustrated in the reference pane of Fig. 9. Wavelet transformations fail to sparsify this wavefield sufficiently because of pervasive high-frequency content [e.g., n(20) = 0.793 for the Daubechies 16 wavelet transformation], which results in reconstructions that degrade rapidly as the undersampling factor decreases.In particular, the tips gradually disappear (see the amplitude images in Fig. 9) and the phase gets corrupted; first at the tips and then towards the bulk as the undersampling factor decreases.In order to avoid this stumbling block, we have to exclude also wavefields whose phase changes too fast from one pixel to another.

Experiments with regular object wavefields
Since the sparsity metric n(p) is also used in thresholding methods for wavelet denoising [47], it is not sensitive to additive noise as long as p is sufficiently small.However, p should also be reasonably large such that the most significant wavelet coefficients represent the actual object wavefield accurately.Table 3 shows the values of n(20) and n(30) for the 3 noisy object wavefields and eleven wavelet transformations.We expect for the USAF resolution chart values of the same order as those in Table 2. Following this assumption we select n(20) as the sparsity metric for the noisy object wavefields.Furthermore, we set the undersampling factor δ to 0.1 and the residual threshold [see Eq.  3 indicates that the Haar wavelet transformation is not the most sparsifying wavelet transformation for the PWC USAF resolution chart.For the orthogonal wavelet transformations, the Daubechies 8 and Battle-Lemarié 3 wavelet transformations are the most sparsifying operators, closely followed by the Daubechies 16 wavelet transformation.We note that these findings contradict with the results of the simulations; probably caused by the Gaussian shape of the laser's TEM00 mode and object features that are smoothed by the optical system.As a consequence, the sparsifying operator needs more vanishing moments.Visual inspection of  In Table 3, we find that the CDF 17/11 and reverse CDF 9/7 wavelet transformations yield the sparsest representations among all wavelet transformations.The values of n(20) are nevertheless very close to each other and comparable with the results for the Daubechies 8, Daubechies 16, and Battle-Lemarié 3 wavelet transformations, just as it is the case for the simulations.For the two most sparsifying CDF wavelet transformations, the amplitude images of the reconstructions are shown in Figs.11(a) and 11(b).we discern that the numbers on the USAF resolution chart and the grains in the leaf are reconstructed with at least the same quality as for any of the orthogonal wavelet transformations in Fig. 10.Our experiments show thus that also these specific biorthogonal wavelet transformations are suited as sparsifying operators in CDH.
Finally, the values of n(20) for the B-spline 3 wavelet transformations are rather high in comparison with the other wavelet transformations, see Table 3.More dense representations in combination with weaker undersampling bounds (see Fig. 2) result in reconstructions that are severly distorted.Another argument that can directly explain these distortions is that the restricted isometry constants for B-spline 3 wavelet transformations are the largest among all wavelet transformations we investigated, resulting in less robustness to noise and sparsity defects.Not all biorthogonal wavelet transformations are thus suited as sparsifying operator in CDH.
As a final note, we draw the attention again to the number of vanishing moments.For the compression of photographic images, the JPEG 2000 encoder uses CDF 5/3 and CDF 9/7 wavelet transformations with respectively 2 and 4 vanishing moments.We recommend for the use of wavelet transformations in holographic imaging methodologies 4 to 8 vanishing moments because the complex-valued nature of wavefields brings more and richer features along.
We note that we have to resort to another solver when dealing with the total variation operator, namely NESTA [48].
The reconstructions are shown in Fig. 12.We observe that the discrete cosine transformation leads to very granular reconstructions: this transformation simply yields a representation that is not sparse enough.For the case of total variation minimization, the results are content dependent.The USAF resolution chart is reconstructed better than with CDF wavelet transformations as sparsifying operator; the bars up to element 2 of group 5 are distinguishable, nicely contoured, and the edges and corners are very sharp.For the leaf, both the reconstructions using CDF wavelet transformations and the total variation minimizer show artefacts.The amplitude of the reconstruction of the former shows small dots whereas the latter is subject to what looks like brush strokes.The phase, however, looks slightly better contoured when using total variation minimization.Finally, the edge above the "15 • " symbol is sharper for the CDF wavelet transformations that for total variation minimization.Depending on the content in the scene, we find thus that the CDF wavelet transformations perform slightly worse or better than a total variation minimization prior.

Fig. 2 .
Fig. 2. Phase transition diagram for CDH with sparsity in a basis of Haar wavelets.All signals are decomposed at 9 scales, d = 1, N = 512, λ = 633 nm, ∆x = 9.3 µm, and z = N (∆x) 2 /λ.The blue and orange isolines are the phase transition bounds for respectively the dual B-spline 3 and B-spline 3 wavelet transformation.See Data File 1 for underlying values.

Fig. 5 .
Fig. 5. Reconstruction error for the computer-generated object wavefields of the USAF resolution chart (b)] and the lens-shaped reflective surface [(c), (d)] as a function of the undersampling factor.The wavelet transformations in Table 2 are deployed as sparsifying operator.See Data File 2 for underlying values.

Fig. 6 .
Fig. 6.Iteration count as a function of the undersampling factor δ for (a) the USAF resolution chart, (b) the lens-shaped relfective surface, (•) the Haar wavelet transformation, and ( ) the Battle-Lemarié 3 wavelet transformation.The reconstruction error is visualized by the color of the markers.

Fig. 9 .
Fig. 9.The reconstructions of the cantilevers degrade from the tips (rapid phase variations) towards the support (slower phase variations) as the undersampling factor decreases.The corresponding holograms were captured at a distance of 4.48 cm from the cantilevers.

( 3 )
] to 0.02 times the Euclidean norm of the undersampled hologram to account for additive noise and other disturbances (increase by 1 % until there is convergence in a few hundred iteration steps).Regarding the reconstruction errors, we prefer visual inspection over the computation of SERs because there are no ground-truth references for the object wavefields.Figs. 10 and 11 depict the amplitude of the reconstructions for respectively 3 orthogonal and biorthogonal wavelet transformations.At first, Table

Fig. 10 (
Fig. 10(b) reveals blocking artifacts for the Haar wavelet transformation whereas the other reconstructions are not affected by this type of distortion, resulting in smoother and more detailed reconstructed object wavefields.In Table3, we find that the CDF 17/11 and reverse CDF 9/7 wavelet transformations yield the sparsest representations among all wavelet transformations.The values of n(20) are nevertheless very close to each other and comparable with the results for the Daubechies 8, Daubechies 16, and Battle-Lemarié 3 wavelet transformations, just as it is the case for the simulations.For the two most sparsifying CDF wavelet transformations, the amplitude images of the reconstructions are shown in Figs.11(a) and 11(b).we discern that the numbers on the USAF resolution chart and the grains in the leaf are reconstructed with at least the same quality as for any of the orthogonal wavelet transformations in Fig.10.Our experiments show thus that also these specific biorthogonal wavelet transformations are suited as sparsifying operators in CDH.Finally, the values of n(20) for the B-spline 3 wavelet transformations are rather high in comparison with the other wavelet transformations, see Table3.More dense representations in combination with weaker undersampling bounds (see Fig.2) result in reconstructions that are severly distorted.Another argument that can directly explain these distortions is that the restricted isometry constants for B-spline 3 wavelet transformations are the largest among all wavelet transformations we investigated, resulting in less robustness to noise and sparsity defects.Not all biorthogonal wavelet transformations are thus suited as sparsifying operator in CDH.As a final note, we draw the attention again to the number of vanishing moments.For the compression of photographic images, the JPEG 2000 encoder uses CDF 5/3 and CDF 9/7 wavelet transformations with respectively 2 and 4 vanishing moments.We recommend for the use of wavelet transformations in holographic imaging methodologies 4 to 8 vanishing moments because the complex-valued nature of wavefields brings more and richer features along.

6. 3 .
Discrete cosine transformation and total variation minimizationIn the final experiment, we compare the reconstructions for the CDF wavelet transformations with the reconstructions using the discrete cosine transformation as sparsifying operator and the solution of the total variation minimization problem α = min α TV(α) subject to g T − TΦΨα 2 < ε.

Fig. 12 .
Fig. 12. Reconstructed wavefields of a USAF resolution chart (amplitude; first column), a leaf (amplitude; second column) (phase; third column), and an angle grid (amplitude; fourth column).The holograms are undersampled 10 times (δ = 0.1).The sparsifying operator for row (a) is the discrete cosine transform and for row (b) total variation minimization is used as prior information.

Table 1 .
Properties of 11 wavelet transformations wavelet transformation short name orthogonality a V columns of A T a ⊥ stands for "orthogonal".b Each wavelet is scaled such that its Euclidean norm is equal to 1; this procedure improves the results significantly.

Table 2 .
Sparsity of the computer-generated object wavefields in eleven wavelet bases a a For every column, the normalized sparsity metrics n that are underlined / typed in bold indicate that the corresponding wavelet transformation is the most sparsifying one among its class members / all listed wavelet transformations. .

Table 3 .
Sparsity of the experimental object wavefields in eleven wavelet domains a For every column, the values that are underlined / typed in bold indicate that the corresponding wavelet transformation is the most sparsifying one among its class members / all listed wavelet transformations. a