Depth-resolved holographic reconstructions by three-dimensional deconvolution

: Methods of three-dimensional deconvolution with a point-spread function as frequently employed in optical microscopy to reconstruct true three-dimensional distribution of objects are extended to holographic reconstructions. Two such schemes have been developed and are discussed: an instant deconvolution using the Wiener filter as well as an iterative deconvolution routine. The instant 3d-deconvolution can be applied to restore the positions of volume-spread objects such as small particles. The iterative deconvolution can be applied to restore the distribution of complex and extended objects. Simulated and experimental examples are presented and demonstrate artifact and noise free three-dimensional reconstructions from a single two-dimensional holographic record. three-dimensional Abstract: Methods of three-dimensional deconvolution with a point-spread function as frequently employed in optical microscopy to reconstruct true three-dimensional distribution of objects are extended to holographic reconstructions. Two such schemes have been developed and are discussed: an instant deconvolution using the Wiener filter as well as an iterative deconvolution routine. The instant 3d-deconvolution can be applied to restore the positions of volume-spread objects such as small particles. The iterative deconvolution can be applied to restore the distribution of complex and extended objects. Simulated and experimental examples are presented and demonstrate artifact and noise free three-dimensional reconstructions from a single two-dimensional holographic record.


Introduction
The beauty of holography lies within the possibility to restore three-dimensional information from a two-dimensional interference pattern. In digital holography [1][2][3], the widely used scientific application of holography, the hologram is usually formed as an interference image captured by a CCD camera and has no physical z-thickness. As a result, the three-dimensional distribution of the object wave-front reconstructed from such a digital hologram suffers from a limited depth of focus and the reconstructed objects do not appear well localized but are prolonged in z-direction. Some discussion on depth-resolution in digital holography can be found in [4][5][6], including recently proposed method of depth-resolution improvement by intensity averaging [6]. In the process of reconstruction of digital holograms, the threedimensional object wave-front is obtained as a set of two-dimensional complex distributions at various distances from the hologram plane, which is analogous to 'optical sectioning' in optical microscopy. In every such two-dimensional section there is signal from the in-focus part of the object but also blurred signal from the out-of-focus part of the object. Consequently, the object wave-front reconstructed from a hologram does not directly represent the original object's distribution. Here we present methods for retrieving the depthresolved original object's distribution from its hologram.

Image formation as convolution with the point-spread function
The problem of out-of focus signal has been studied in optical microscopy and has a number of solutions by applying the so-called deconvolution methods; overviews of the methods can be found in [7][8][9][10][11]. When images of a three-dimensional object () Or are acquired at several focal planes (optical sections), a three-dimensional volume of measured data () Mr is created: (1) The image of a point scatterer through the same optical system forms the point-spread function (PSF) of the system: In optical microscopy, the PSF is usually recorded experimentally by using a small object, e.g. a microsphere.
If the object can be represented as a sum of point scatterers (condition of linearity): and the optical system images the individual scatterers in the same way (condition of shift invariance), the resulting image can be described as the sum of images of the individual scatterers [12]: Equation (4) can be viewed as a convolution: The deconvolution of the measured three-dimensional data volume () Mr with the known PSF of the system makes it possible to obtain the three-dimensional object () Or .

Formation and reconstruction of a hologram
An important difference between imaging in a conventional microscope and holography is that coherence and interference are the key elements for the latter. The interference between waves scattered by different parts of the object must be accounted for. The wave-front reconstructed from a hologram carries amplitude and phase information and exhibits a complex distribution. The requirement of linearity is fulfilled for the complex optical fields: the total complex optical field created by a set of scatterers is the sum of the complex optical fields created by each scatterer. An additional requirement in the case of holography is that the effects of multiple scattering between different parts of the object must be negligible. To fulfill the requirement of shift invariance we consider plane waves for simplicity (in case of a spherical wave-front geometrical scaling must be applied in addition). When an incident plane wave is scattered by some three-dimensional object () Or , neglecting multiple scattering effects, the resultant object wave-front at any point r can be represented as a convolution of the object function with the impulse response of free space propagation: where the integration is performed over the object's coordinates O r .  , the distribution in the hologram plane is given by: The mathematical routine for the reconstruction of holograms is independent of the nature of the recorded objects and consists of the following two steps. 1. Illumination with the reference wave and 2. Backward propagation to the object-position. In a first step, the multiplication of the hologram with the simulated complex reference wave RS () Ur, provides the fully restored two-dimensional complex distribution of the object wave in the hologram plane: In a second step, the backward propagation of the complex wavefront from the hologram plane to the object's location is calculated by applying the Huygens-Fresnel principle: (8) and the result of the integration provides the reconstructed three-dimensional wave-front of the object wave. In analogy to optical microscopy, the wave-front reconstructed from the hologram of a point-scatterer shall serve as a basis for the PSF: Equations (8) and (9) are usually referred to in literature as reconstruction by convolution with the impulse response of free space propagation (see for instance [13]). Next, we would like to demonstrate, how the object distribution can be extracted from the reconstructed object wavefront by employing a three-dimensional deconvolution with the point-spread function. For the three-dimensional deconvolution procedure there are two three-dimensional distributions of the reconstructed wave-fronts available: O () Ur and P () Ur.

Three-dimensional deconvolution of the reconstructed optical fields
In direct analogy to optical microscope data, the complex optical field reconstructed from a hologram of the object O () Ur can be considered as the measured signal () Mr and the complex optical field reconstructed from a hologram of a point scatterer as the PSF( ) r . We thus obtain the convolution equation, described by Eq. (5), as: where () Or is the unknown three-dimensional object distribution. The simplest way to solve Eq. (10) is to use the convolution theorem (FT here refers to a three-dimensional Fourier transform): However, such a direct approach will not result in the object's three-dimensional distribution for the following reason. The reconstructed optical fields given by Eqs. (8) and (9) can be represented as a three-dimensional convolution with the impulse response of the free space propagation and the convolution theorem can be applied: and: where the three-dimensional Fourier transform is given by: (15) or, in Cartesian coordinates: (16) where ( , , ) q     are coordinates in the Fourier domain. Next, substitution of Eq. (13) and Eq. (14) into Eq. (12) results in elimination of the threedimensional Fourier transforms of the free-space propagators, and only a fraction of the threedimensional Fourier transforms of the two-dimensional images remains: which after integration over the z-coordinate are in detail By substituting the last two equations into Eq. (17), we obtain: Here, obviously, the fraction in the brackets does not have a dependency on  and thus, the backward three-dimensional Fourier transform will always result in a  -function at 0 z  as a result of integration over d . As a consequence, the object function () Or will always be just a two-dimensional distribution in the plane 0 z  . This implies that such direct 3d-deconvolution applied to the reconstructed complex fields would always result in losing the z -information.
We suggest two ways to solve the problem. 1. A 3d-deconvolution of the intensities of the reconstructed fields. 2. An iterative 3d-deconvolution.

Instant 3d-deconvolution of the reconstructed intensities by the Wiener filter
One way to perform an instant 3d-deconvolution is to apply it to the reconstructed intensity distributions instead of the complex amplitudes. Mathematically, when intensities are considered instead of complex amplitudes, it allows avoiding elimination of z and ηdependent terms while performing the division given by Eq. (12). In fact, all deconvolution methods employed in optical microscopy are designed to be applied to the intensity distributions [7][8][9][10]. We consider the reconstructed intensity of the object wave r . In the case of intensities, the requirement of linearity is only fulfilled when the effects of interference between waves from individual point scatterers constituting the object are negligible. This condition, however, is not unrealistic, but very common for instance in the holographic study of particle fields where a large amount of identical particles is spread in a volume [14,15]. The instant 3d-deconvolution by the Wiener filter is calculated as follows: where β is a small addend to avoid division by zero;  is usually much smaller than 2 P () Ur and selected to achieve a good signal-to-noise reconstruction. The result of the 3d-deconvolution is a set of sharp localizations (ideally,  -functions) of individual scatterers constituting the object.

Iterative deconvolution of the reconstructed complex fields
The method of iterative 3d-deconvolution solves Eq. (11) without direct division. The following example shall be considered to illustrate that. Two point scatterers are given: point scatterer A shall be positioned at A (0,0,0) r  and scatterer B at B . The scattered complex wave-fronts are described by: The two holograms of the scatterers are recorded at SD zz  . The complex field reconstructed from the hologram of scatterer A is employed as the PSF of the optical system. The complex field reconstructed from the hologram of scatterer B is the reconstructed object wave-front. When Eq. (21) and Eq. (22) are substituted into Eq. (18), we obtain the following threedimensional distribution in the Fourier domain: Obviously, when substituted in Eq. (12), the direct division of Eq. (24) with Eq. (23) will lead to the elimination of the  -dependency and consequently the loss of z-information. When Eqs. (24) and (23) are substituted into the convolution Eq. (11), by taking into account the three-dimensional Fourier transform of the impulse response of free space propagator: after some simplifications, we arrive at Eq. (11) as: From this, the exact position of scatterer B is obtained by a three-dimensional backward Fourier transformation. A direct division could thus be avoided and Eq. (11) is solved by carefully matching the left and right sides which reveals the recovery of the object distribution. In the next section we verify this theoretical hypothesis by applying it to simulated and experimental examples. Since this iterative deconvolution is applied to the complex optical fields rather than to the intensities, it can also be applied to continuous extended objects where the interference between different parts of the object cannot be neglected anymore.
The iterative loop is adapted from optical microscopy and includes the following steps [16]: where  is a small addend to avoid division by zero. The selection criteria for  have to fulfill the following conditions:  has to be much smaller than the maximal value of Ur , respectively, while P () Ur remains unaltered. The three-dimensional distributions obtained at each iterative step can be multiplied with the apodization function to avoid accumulation of artifacts during repeated three-dimensional Fourier-transformations. The apodization function selected here consists of a three-dimensional spherical filter of radius R, which leaves all pixel values inside the sphere unchanged and sets the ones outside the sphere to zero.

Numerical simulation and reconstruction of holograms
Before applying the routine to simulated and experimental examples, we would like to describe in some detail how the numerical optical field propagation used in the simulation and reconstruction routines is carried out. For the calculation of the complex field propagation we used the well-known algorithm for plane wave propagation (see for example [17]) that consists of the following three steps: (a) Two-dimensional forward Fourier-transform of the complex field at the plane i z .
(b) Multiplication with the simulated complex transfer function: (c) Two-dimensional backward Fourier-transform of the product (a) and (b).
The reconstruction procedure consists of the same steps (a)-(c), but the complex transfer function is the complex conjugated one due to the back propagation: zz  . A PSF was created as the field distribution reconstructed from a hologram of a point scatterer. The hologram of a point scatterer could be either experimentally recorded or simulated numerically. The hologram for the PSF must be of the same size as the hologram of the object and it should be reconstructed for the same z-span. The result of 3d-deconvolution shows the relative shift of the object's parts to the position of the point scatterer.
All images were sampled with 200x200 pixels size. The reconstructions were done at a set of z-distances for 200 steps in total. The resulting three-dimensional data volume exhibited a size of 200x200x200pixels limited by the ability of the employed MATLAB software to operate with arrays of larger size. To calculate the three-dimensional FFT we used the MATLAB function "fftn", which allows FFT of any dimension, together with "fftshift" and "ifftshift" commands to keep the lower frequencies in the center of the complex images.

Three-dimensional deconvolution of simulated holograms
To simulate a three-dimensional object, we used a set of two-dimensional transmission functions at various z-distances. The complex optical field was consequently propagated between each pair of xy-planes. Thus, multiple scattering effects were taken into account to generate a simulated hologram as realistic as possible to check whether multiple scattering effects have an influence on the 3d-deconvolution.

Instant 3d-deconvolution of the reconstructed intensities by the Wiener filter
For the instant 3d-deconvolution the condition of negligible interference between wave-fields scattered from different parts constituting the object must be satisfied. We thus simulated a particle distribution in the shape of the three Greek letters α, β and γ, each made up of point scatterers and placed at three planes at different distances from the screen; see Fig. 2(a).
The letters α, β and γ were positioned at a 90, 80 and 70mm distance from the screen plane, respectively. The incident radiation exhibited a wavelength of 500nm. Object and the screen area sizes were selected to be 2 mm. The reconstructions were carried out from z = 40 to 120mm distance with a step-width of 0.4mm for a number of 200 steps in total. The PSF was created as the reconstruction of the hologram of a point scatterer positioned at 80mm from the screen; see Fig. 2 The three-dimensional amplitude distribution reconstructed from the simulated hologram shows the presence of the in-focus and out-of-focus signal, see Fig. 3. The instant 3d-deconvolution was performed using Eq. (20) with β = 1. The 3d-deconvolution leads to the removal of the out-of-focus signal. The individual scatterers revealed sharply and the patterns of the three Greek letters could clearly be distinguished, see Fig. 3 and Fig. 4. The z-localization of a selected particle in the deconvoluted distribution amounts to less than 1.6mm in the reconstruction volume distributed over 80mm in z-direction, see Fig. 4(d). In addition, any noise caused by the out-of focus signals and twin images was removed.

Iterative 3d-deconvolution of the reconstructed complex field
The iterative 3d-deconvolution was applied to the same simulated hologram of the particle distribution using Eqs. Ur . The result after 24 iterations is shown in Fig. 5. The z-localization of a selected particle amounts to just 1.6mm in the deconvoluted object distribution, as shown in Fig. 5(b). The error function, computed as: shows a decrease after the first iterations followed by oscillations around some minimal value. However, the error function does not actually reach zero. This is because the updated complex field (k) O () Ur which is simulated as the convolution of the updated object distribution function (k) () Or with the PSF does not include such effects as multiple scattering or twin images and is thus never exactly equal to the reconstructed complex field O () Ur .   The iterative 3d-deconvolution can also be applied to restore continuous, non-point-like objects. As an example, we simulated a hologram of the three continuous Greek letters: α, β and γ placed at a 90, 80 and 70mm distance from the screen. The other parameters of the hologram simulation and reconstruction remained the same as in the example described above. The 3d-deconvolution was done with the iterative loop Eqs. (28)-(29) with β = 0.1. The function (k 1) () Or  obtained after every step (ii) in the iterative loop was multiplied with the apodization function -a three-dimensional spherical filter of radius R = 80 pixels. To smooth the appearance of continuous objects and avoid divergence of the routine caused by random peaks of noise, a Gaussian low-pass filter was applied: where μ,ν,η are coordinates in the Fourier-domain and D is a small constant. In our example, the filter was applied in every 5th iteration with D = 5. The results of the iterative 3d-deconvolution are shown in Fig. 6.

Experimental optical holograms
To apply our method also to experimental data, we recorded optical holograms of a threedimensional object: a piece of a 170μm thin glass plate with polystyrene microspheres of 4μm in diameter randomly distributed on both sides of the glass (S1 and S2, see Fig. 7). The employed holographic inline scheme using plane waves is shown in Fig. 7. The wavelength of the green laser employed amounts to 532nm, the distance between object and microscope objective was about 1mm, the imaged area about 550μm in diameter, the distance between the microscope objective and the screen was selected to 16cm, a hologram of about 30mm in diameter was recorded on a semitransparent screen and captured by a CCD camera. Fig. 7. Experimental scheme for in inline holographic recording with plane waves.
Directly after the acquisition of the hologram H0 of the sample, an image in the absence of the sample was recordedthe background image B. Thereafter the laser beam was blocked by a piece of black cardboard and the intensity of the dark environment -IDEwas recorded.
The IDE amounts to about 50 gray levels on the 10 Bit (1024 gray levels) CCD camera chip. The normalization of the hologram was calculated by the division H = (H0-IDE)/(B-IDE). The result was a contrast hologram H of average intensity of 1; an example of which is shown in Fig. 8. Polystyrene exhibits an absorption coefficient of about 3cm 1 [18] and a refractive index of about 1.6 at the wavelength of 532nm. A 4μm diameter polystyrene sphere thus constitutes a weak absorbing but a strong (up to 28 radians) phase object. Below we show that even such strong phase shifting properties represent no problem for our deconvolution routine.
The recorded holograms exhibited an original size of 1000x1000 pixels but needed to be re-sampled to 200x200 pixels for the reconstruction and 3d-deconvolution process for reasons explained in Section 5. Re-sampling was done using a linear interpolation in both, x-and ydirection, for computing pixel locations and their gray values. The information about individual scatterers is stored in the interference fringes which are distributed over the entire hologram surface. Thus, re-sampling of the hologram does not alter the information about the position of the scatterers but just introduces a loss of spatial resolution in the reconstructions.

Instant 3d-deconvolution of the reconstructed intensities using the Wiener filter
The optical holograms of the spheres were normalized and reconstructed, see Fig. 9. The reconstructions were carried out from z = 170 to 1770μm distance with a step-width of 8μm for a number of 200 steps in total. The two planes of in-focus spheres were found at distances of 970μm (S1) and 800μm (S2) from the hologram plane, as shown in Fig. 9(d). The hologram for the PSF was created as a cut-out from a hologram of the spheres where just an individual sphere was detected, see Fig. 9(a), and the cut-out was subsequently placed into the image's center and zero-padded, as shown in Fig. 9(b). The instant 3d-deconvolution using Eq. (20) with β = 0.01 was applied to the reconstructed intensity distribution. As a result, the well-localized z-positions of the objects were retrieved free from out-of focus-signals, as shown in Fig. 9(e).

Iterative 3d-deconvolution of the reconstructed complex fields
For the iterative 3d-deconvolution we used the complex distribution of the optical field reconstructed from the same experimental hologram as in the example above. The iterative 3d-deconvolution was done using 12 iterations by applying Eqs. (28)-(29) with β = 0.012. The function ( 1) obtained after every step (ii) in the iterative loop was multiplied with the apodization function -a three-dimensional spherical filter of radius R = 98 pixels. No other numerical filters were applied. The results are shown in Fig. 9(f) and exhibit superior removal of the out-of-focus signal when compared to that of the instant 3d-deconvolution using the Wiener filter, as evident from comparing Fig. 9(e) with Fig. 9(f).
Apparently, the 3d-deconvolution also improves the resolution in the xy-plane, as illustrated in Fig. 10 where we show the magnified images of the reconstructions at the plane S1 before and after the iterative 3d-deconvolution. While the clusters of spheres are originally not well resolved, after applying the 3d-deconvolution scheme, the positions of the individual spheres composing the clusters are clearly visible and emerge as individual sharp spots.

Three-dimensional deconvolution of the reconstructed optical fields with a simulated PSF
In some cases, the PSF is not available from experimental data. For instance, when a hologram does not contain an image of an individual stand-alone sphere. In this case the hologram of a point scatterer can be simulated and its complex reconstruction provides the PSF. To check this approach, we simulated a hologram of a point scatterer in the middle of a 550x550μm 2 area placed at 970μm distance from the screen, as shown in Fig. 11(a) and we used the reconstruction of its hologram as the PSF for 3d-deconvolution. For the instant 3d-deconvolution we selected β = 0.01 and for the iterative 3d-deconvolution β = 0.2. The results of the instant and iterative 3d-deconvolutions are shown in Fig. 11(b) and 11(c), respectively. The resulting distributions of spheres are similar to the distributions obtained by 3d-deconvolution with the experimental PSF, shown in Fig. 9(e) and 9(f). The instant 3d-deconvolution using the simulated PSF, as shown in Fig. 11(b), results in point-like reconstructed objects, however with some remaining background noise left. The iterative 3d-deconvolution with the simulated PSF provides noise-free reconstructed spheres, however they appear slightly extended in z-direction. The result demonstrates that comparably good results of 3d-deconvolution can be obtained using either an experimentally obtained PSF or a simulated one.    Reconstructed amplitude distribution at planes S1 and S2 and its three-dimensional representation. (c) Results of the iterative 3d-deconvolution at the plane S1 and S2 together with its three-dimensional representations. (d) Amplitude profiles along z-direction at the spheres' positions before (black) and after (red) 3d-deconvolution.

Iterative 3d-deconvolution of the reconstructed complex fields applied to extended objects
In this section we show that the method of the iterative 3d-deconvolution with the simulated PSF can also be applied to complex and/or not-point-like objects, such as two spheres of a extended size. A 200x200pixel area corresponding to 110x110μm 2 was selected in an experimental hologram of spheres; see Fig. 12(a). The reconstructions were carried out from z = 10 to 1610μm distance with a step-width of 8μm for a number of 200 steps in total. The amplitude reconstruction shows two spheres at different z-distances: at 810μm (S1) and 630μm (S2) from the hologram. The PSF was obtained by the reconstruction of the complex field distribution from the simulated hologram of a point scatterer placed in the middle of 110x110μm 2 area at the distance of z = 810μm from the screen. The iterative 3ddeconvolution was performed with 12 iterations using β = 0.012 without any numerical filters. The results are shown in Figs. 12(b)-12(d): after the iterative 3d-deconvolution the shape of the spheres appears sharper in all three directions and the out-of focus signal is removed.

Conclusions
We have shown that 3d-deconvolution methods as known from optical microscopy can also be applied to restore true object distribution from holographic records and their reconstructions. We demonstrated methods of instant and iterative 3d-deconvolutions applied to simulated as well as to experimental holograms. The three-dimensional distribution of polystyrene spheres deposited on glass surfaces could be retrieved free from disturbing out-of focus and twin image signals by a 3d-deconvolution of the hologram reconstructions. Instant 3d-deconvolution can be applied if the object consists of a set of spaced-apart particles distributed in a certain volume. The method of instant 3d-deconvolution has potential to be applied in holographic particle tracking velocimetry, where the identification of the correct z-position of the object and the removal of the out-of-focus signal is essential. The iterative 3d-deconvolution can also be applied for objects exhibiting complex and/or extended shapes or distributions. In view of future applications, it is important to note that the methods of 3d-deconvolution are independent of the nature of the coherent radiation used. This opens up the possibility for 3d reconstruction also in X-ray or electron holography. The PSF required for 3d-deconvolutions can either be obtained by experimentally recording and reconstructing a hologram of an arbitrarily small individual scatterer or as the reconstruction of a simulated hologram of a point scatterer, which in general is similar to a Fresnel zoneplate. The approach of creating the PSF from a simulated hologram can also be applied to account for angular-dependent scattering, as it is the case of low-energy electron scattering. A realistic but simulated PSF can thus be used in electron holography where an experimental hologram of a point scatterer is not readily available.
By employing 3d-deconvolution methods the out-of focus signal is brought back to its scatterer and the twin images are automatically removed as they are not part of the scattered wave. Thus, spatially well localized parts of an object are recovered free from artifacts and the beauty of holography as a three-dimensional imaging technique lives up to its full potential.