Depth-resolved holographic reconstructions by three-dimensional deconvolution

Tatiana Latychevskaia; Fabian Gehri; Hans-Werner Fink

doi:10.1364/OE.18.022527

1. Introduction

The beauty of holography lies within the possibility to restore three-dimensional information from a two-dimensional interference pattern. In digital holography [1–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–6], including recently proposed method of depth-resolution improvement by intensity averaging [6]. In the process of reconstruction of digital holograms, the three-dimensional 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 depth-resolved original object’s distribution from its hologram.

2. 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–11]. When images of a three-dimensional object $O (\vec{r})$ are acquired at several focal planes (optical sections), a three-dimensional volume of measured data $M (\vec{r})$ is created:

O (\vec{r}) \Rightarrow M (\vec{r}) .

The image of a point scatterer through the same optical system forms the point-spread function (PSF) of the system:

δ (\vec{r}) \Rightarrow PSF (\vec{r}) .

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):

O (\vec{r}) = \int O (\vec{s}) δ (\vec{r} - \vec{s}) d \vec{s}

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]:

M (\vec{r}) = \int O (\vec{s}) PSF (\vec{r} - \vec{s}) d \vec{s} .

Equation (4) can be viewed as a convolution:

M (\vec{r}) = O (\vec{r}) \otimes PSF (\vec{r}) .

The deconvolution of the measured three-dimensional data volume

M (\vec{r})

with the known PSF of the system makes it possible to obtain the three-dimensional object

O (\vec{r})

.

3. 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 $O (\vec{r})$ , neglecting multiple scattering effects, the resultant object wave-front at any point $\vec{r}$ can be represented as a convolution of the object function with the impulse response of free space propagation:

U_{Object} (\vec{r}) = ∭ O ({\vec{r}}_{O}) \frac{\exp (i k | {\vec{r}}_{O} - \vec{r} |)}{| {\vec{r}}_{O} - \vec{r} |} d {\vec{r}}_{O} = O (\vec{r}) \otimes \frac{\exp (i k r)}{r}

where the integration is performed over the object’s coordinates

{\vec{r}}_{O}

.

The superposition of the wave scattered by the object $U_{Object} ({\vec{r}}_{S})$ and a reference wave $U_{R} ({\vec{r}}_{S})$ , where ${\vec{r}}_{S} = (x_{S}, y_{S}, z_{S})$ , when recorded at some distance $z_{D}$ produces a hologram, see Fig. 1 . The distribution of the transmission function of the recorded hologram is given by $| \tilde{H} ({\vec{r}}_{S}) | = {| {\tilde{U}}_{Object} ({\vec{r}}_{S}) + {\tilde{U}}_{R} ({\vec{r}}_{S}) |}^{2}$ , where the tilde denotes the two-dimensionality, i.e. the function is defined only in the screen $(x_{S}, y_{S})$ -plane. After subtraction of the constant background ${| {\tilde{U}}_{R} ({\vec{r}}_{S}) |}^{2}$ and in the approximation of a weak object wave ${| {\tilde{U}}_{Object} ({\vec{r}}_{S}) |}^{2} < < {| {\tilde{U}}_{R} ({\vec{r}}_{S}) |}^{2}$ , the distribution in the hologram plane is given by:

\tilde{H} ({\vec{r}}_{S}) ~ {\tilde{U}}_{Object} ({\vec{r}}_{S}) {\tilde{U}}_{R}^{*} ({\vec{r}}_{S}) + {\tilde{U}}_{Objext}^{*} ({\vec{r}}_{S}) {\tilde{U}}_{R} ({\vec{r}}_{S}) .

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

{\tilde{U}}_{R} ({\vec{r}}_{S})

, provides the fully restored two-dimensional complex distribution of the object wave in the hologram plane:

{\tilde{U}}_{R} ({\vec{r}}_{S}) \tilde{H} ({\vec{r}}_{S}) ~ {\tilde{U}}_{Object} ({\vec{r}}_{S})

. In a second step, the backward propagation of the complex wave-front from the hologram plane to the object's location is calculated by applying the Huygens-Fresnel principle:

U_{O} (\vec{r}) = ∭ {\tilde{U}}_{Object} ({\vec{r}}_{S}) δ (z_{S} - z_{D}) \frac{\exp (i k | {\vec{r}}_{S} - \vec{r} |)}{| {\vec{r}}_{S} - \vec{r} |} d {\vec{r}}_{S}

and the result of the integration provides the reconstructed three-dimensional wave-front of the object wave.

Fig. 1 Scheme of hologram recording and illustration to the used symbols.

Download Full Size | PDF

In analogy to optical microscopy, the wave-front reconstructed from the hologram of a point-scatterer shall serve as a basis for the PSF:

U_{P} (\vec{r}) = ∭ {\tilde{U}}_{Point} ({\vec{r}}_{S}) δ (z_{S} - z_{D}) \frac{\exp (i k | {\vec{r}}_{S} - \vec{r} |)}{| {\vec{r}}_{S} - \vec{r} |} d {\vec{r}}_{S} .

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 wave-front 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:

U_{O} (\vec{r})

and

U_{P} (\vec{r})

.

4. 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 $U_{O} (\vec{r})$ can be considered as the measured signal $M (\vec{r})$ and the complex optical field reconstructed from a hologram of a point scatterer as the $PSF (\vec{r})$ . We thus obtain the convolution equation, described by Eq. (5), as:

U_{O} (\vec{r}) = O (\vec{r}) \otimes U_{P} (\vec{r}),

where

O (\vec{r})

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):

FT(U_{O} (\vec{r})) = FT (O (\vec{r})) FT(U_{P} (\vec{r}))
.

and solve the equation by division:

O (\vec{r}) = {FT}^{- 1} (\frac{FT(U_{O} (\vec{r}))}{FT(U_{P} (\vec{r}))}) .

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:

FT(U_{O} (\vec{r}))=FT({\tilde{U}}_{Object} (\vec{r}) δ (z - z_{D}))FT (\frac{\exp (i k r)}{r})

and:

FT(U_{P} (\vec{r}))=FT({\tilde{U}}_{Point} (\vec{r}) δ (z - z_{D}))FT (\frac{\exp (i k r)}{r}),

where the three-dimensional Fourier transform is given by:

FT (U (\vec{r})) = ∭ U (\vec{r}) \exp (- 2 π i \vec{r} \vec{q}) d \vec{r}

or, in Cartesian coordinates:

FT (U (x, y, z)) = ∭ U (x, y, z) \exp (- 2 π i (x μ + y ν + z η)) d x d y d z

where

\vec{q} = (μ, ν, η)

are coordinates in the Fourier domain.

Next, substitution of Eq. (13) and Eq. (14) into Eq. (12) results in elimination of the three-dimensional Fourier transforms of the free-space propagators, and only a fraction of the three-dimensional Fourier transforms of the two-dimensional images remains:

O (\vec{r}) = {FT}^{- 1} (\frac{FT({\tilde{U}}_{Object} (\vec{r}) δ (z - z_{D}))}{FT({\tilde{U}}_{Point} (\vec{r}) δ (z - z_{D}))}),

which after integration over the z-coordinate are in detail

FT({\tilde{U}}_{Object} (\vec{r}) δ (z - z_{D})) = exp(-2 π i z_{D} η) \iint {\tilde{U}}_{Object} (x, y)exp(-2 π i (x μ + y ν))d x d y

By substituting the last two equations into Eq. (17), we obtain:

O (\vec{r}) = {FT}^{- 1} (\frac{\iint {\tilde{U}}_{Object} (x, y)exp(-2 π i (x μ + y ν))d x d y}{\iint {\tilde{U}}_{Point} (x, y)exp(-2 π i (x μ + y ν))d x d y}) .

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

z = 0

as a result of integration over

d η

. As a consequence, the object function

O (\vec{r})

will always be just a two-dimensional distribution in the plane

z = 0

. 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.

4.1 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–10]. We consider the reconstructed intensity of the object wave ${| U_{O} (\vec{r}) |}^{2}$ as the measured signal $M (\vec{r})$ and the reconstructed intensity of a point scatterer ${| U_{P} (\vec{r}) |}^{2}$ as the $PSF (\vec{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:

O (\vec{r}) = {FT}^{- 1} (\frac{FT ({| U_{O} (\vec{r}) |}^{2})}{FT ({| U_{P} (\vec{r}) |}^{2} + β)}),

where β is a small addend to avoid division by zero; β is usually much smaller than

{| U_{P} (\vec{r}) |}^{2}

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.

4.2 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 ${\vec{r}}_{A} = (0, 0, 0)$ and scatterer B at ${\vec{r}}_{B} = (x_{B}, y_{B}, z_{B})$ . The scattered complex wave-fronts are described by:

U_{A} ({\vec{r}}_{S}) \approx \frac{\exp (i k \sqrt{x_{S}^{2} + y_{S}^{2} + z_{S}^{2}})}{z_{S}}

and

U_{B} ({\vec{r}}_{S}) \approx \frac{\exp (i k \sqrt{{(x_{S} - x_{B})}^{2} + {(y_{S} - y_{B})}^{2} + {(z_{S} - z_{B})}^{2}})}{(z_{S} - z_{B})} .

The two holograms of the scatterers are recorded at

z_{S} = z_{D}

. 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 three-dimensional distribution in the Fourier domain:

FT({\tilde{U}}_{A} (\vec{r}) δ (z - z_{D})) = exp(-2 π i z_{D} η) \exp (- π i λ z_{D} (μ^{2} + ν^{2}))

FT({\tilde{U}}_{B} (\vec{r}) δ (z - z_{D})) =

= exp(-2 π i z_{D} η)exp(-2 π i (x_{B} μ + y_{B} ν))exp(- π i λ (z_{D} - z_{B}) (μ^{2} + ν^{2}))
.

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:

FT (\frac{\exp (i k r)}{r}) \approx δ (η + \frac{λ}{2} (μ^{2} + ν^{2})),

after some simplifications, we arrive at Eq. (11) as:

exp(-2 π i (x_{B} μ + y_{B} ν))exp(π i λ z_{B} (μ^{2} + ν^{2})) δ (η + \frac{λ}{2} (μ^{2} + ν^{2})) =

= FT (O (\vec{r})) δ (η + \frac{λ}{2} (μ^{2} + ν^{2})) .

(26)

η = - \frac{λ}{2} (μ^{2} + ν^{2})

FT (O (\vec{r})) = \exp (- 2 π i (x_{B} μ + y_{B} ν + z_{B} η)) .

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]:

O^{(1)} (\vec{r}) = U_{O} (\vec{r})

U_{O}^{(k)} (\vec{r}) = O^{(k)} (\vec{r}) \otimes U_{P} (\vec{r})

O^{(k + 1)} (\vec{r}) = O^{(k)} (\vec{r}) \frac{U_{O} (\vec{r})}{U_{O}^{(k)} (\vec{r})}

k = k + 1

The division in step (ii) is calculated as

O^{(k + 1)} (\vec{r}) = O^{(k)} (\vec{r}) \frac{U_{O} (\vec{r}) {(U_{O}^{(k)} (\vec{r}))}^{*}}{{| U_{O}^{(k)} (\vec{r}) |}^{2} + β}

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

{| U_{O}^{(k)} (\vec{r}) |}^{2}

but sufficiently large to ensure convergence of the iterative process. With each iteration step, the object function

O^{(k)} (\vec{r})

and the reconstructed object wave

U^{(k)}_{O} (\vec{r})

are approaching the complex distribution of the true object function

O (\vec{r})

and the reconstructed wave-front

U_{O} (\vec{r})

, respectively, while

U_{P} (\vec{r})

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.

5. 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 $z_{i}$ .
(b) Multiplication with the simulated complex transfer function: $\exp (\frac{2 π i}{λ} z \sqrt{1 - μ^{2} - ν^{2}})$ where $z = z_{i+1} - z_{i}$ is the propagation distance, $μ = \frac{λ}{S} p$ , $ν = \frac{λ}{S} q$ , where S is the area size, $p= - \frac{N}{2} ... \frac{N}{2}$ and $q= - \frac{N}{2} ... \frac{N}{2}$ are the pixel indexes, and N is the number of pixels.
(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:

z \to - z

.

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.

6. 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.

6.1 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) .

Fig. 2 Arrangement of numerical experiments for the simulation of holograms of a particle distribution and of an individual point scatterer. (a) Simulation of the hologram of a particle distribution showing the transmission functions in the three planes at which the particles are arranged. (b) Simulation of the hologram of a point scatterer.

Download Full Size | PDF

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(b).

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.

Fig. 3 (a) Three-dimensional representation of the amplitude of the reconstructed complex wave-front. (b) Amplitude of the complex PSF. (c) Result of the instant 3d-deconvolution.

Download Full Size | PDF

Fig. 4 Results of the instant 3d-deconvolution of the intensities reconstructed from the simulated hologram of a particle distribution. (a) Original objects used in the simulation at the three z-planes. (b) Reconstructed amplitude of the object wave-front at the same z-planes. (c) Result of instant 3d-deconvolution. (d) Intensity profiles along z-direction for the particle indicated by the red arrow in (c) before (black) and after (red) instant 3d-deconvolution.

Download Full Size | PDF

6.2 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. (28)-(29) with β = 0.01. After each step (i) of the iterative loop, the absolute value of the function $U_{O}^{(k)} (\vec{r})$ was normalized so that its minimum and maximum match the absolute value of $U_{O} (\vec{r})$ . 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:

Error (k) = \frac{\sum_{i,j,l} | | U_{O}^{(k)} (i,j,l) | - | U_{O} (i,j,l) | |}{\sum_{i,j,l} | U_{O} (i,j,l) |}

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

U^{(k)}_{O} (\vec{r})

which is simulated as the convolution of the updated object distribution function

O^{(k)} (\vec{r})

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

U_{O} (\vec{r})

.

Fig. 5 Results of the iterative 3d-deconvolution of the complex fields reconstructed from the simulated hologram of a particle distribution. (a) Amplitude distributions of the deconvoluted fields at the three planes. (b) Intensity profiles for a 16mm span around the in-focus position along z-direction for the particle indicated by an arrow in (a) before and after deconvolution. (c) Error as a function of the iteration number.

Download Full Size | PDF

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 $O^{(k + 1)} (\vec{r})$ 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:

GLPF (μ, ν, η) = \exp (- \frac{μ^{2} + ν^{2} + η^{2}}{2 D^{2}})

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 .

Fig. 6 Results of the iterative 3d-deconvolution of the complex fields reconstructed from the simulated hologram of continuous objects. (a) Original objects used in the simulation at the three different planes. (b) Reconstructed objects at the same planes. (c) Result of the 3d-deconvolution. (d) Three-dimensional representation of the reconstructed amplitude of the object wave-front and amplitude of the deconvoluted object function. (e) Intensity profiles for a 16mm span around the in-focus position along z-direction for the point of the γ-object indicated by the red arrow in (c) before (black) and after (red) iterative 3d-deconvolution.

Download Full Size | PDF

7. Experimental optical holograms

To apply our method also to experimental data, we recorded optical holograms of a three-dimensional 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.

Download Full Size | PDF

Directly after the acquisition of the hologram H0 of the sample, an image in the absence of the sample was recorded – the background image B. Thereafter the laser beam was blocked by a piece of black cardboard and the intensity of the dark environment – IDE – was 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 .

Fig. 8 Experimentally recorded hologram H0 of polystyrene spheres (left), the background image B (center) and the normalized hologram H (right).

Download Full Size | PDF

Polystyrene exhibits an absorption coefficient of about 3cm⁻¹ [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 y-direction, 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.

7.1 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).

Fig. 9 3d-deconvolution of experimental holographic reconstruction. (a) Normalized hologram. (b) Hologram for the PSF as a cut out of the experimental hologram related to just a single sphere. (c) Amplitude reconstruction of the hologram for the PSF, slice in xz-plane. (d) Amplitude reconstruction at both sides of the microscopic slide, at the S1 and S2 plane. In-focus and out-of of focus spheres can be distinguished. Right: 3d-representation of the reconstructed amplitude distribution. (e) Results after instant 3d-deconvolution. Amplitude distributions at the S1 and S2 plane together with the 3d-distribution. (f) Results after iterative 3d-deconvolution. Amplitude distributions at the S1 and S2 plane and 3d-distribution.

Download Full Size | PDF

7.2 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 $O^{(k + 1)} (\vec{r})$ 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.

Fig. 10 Magnified images of the reconstructed amplitude before and after iterative 3d-deconvolution at plane S1.

Download Full Size | PDF

7.3 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² 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.

Fig. 11 3d-deconvolution of the experimental holographic reconstruction with a simulated PSF. (a) Simulated hologram of a point scatterer. (b) Results of instant 3d-deconvolution at S1 and S2 planes and as 3d-representation. (b) Results of iterative 3d-deconvolution at S1 and S2 planes and as 3d-representation.

Download Full Size | PDF

7.4 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² 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² area at the distance of z = 810μm from the screen. The iterative 3d-deconvolution 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.

Fig. 12 Results of iterative 3d-deconvolution of the optical field reconstructed from an experimental hologram with a simulated PSF. (a) Selected fragment of an experimental hologram of spheres. (b) 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.

Download Full Size | PDF

8. 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 zone-plate. 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.

Acknowledgements

The authors are grateful to Conrad Escher for sample preparation and to Daniel Eaton for careful reading of the manuscript. Financial support of the European Project SIBMAR and the Forschungskredit of the University of Zurich are gratefully acknowledged.

References and links

1. J. W. Goodman and R. W. Lawrence, “Digital image formation from electronically detected holograms,” Appl. Phys. Lett. 11(3), 77–79 (1967). [CrossRef]

2. L. P. Yaroslavsky, and N. S. Merzlykaov, Methods of digital holography, (Consultants Bureau, New York, 1980).

3. U. Schnars, and W. Jueptner, Digital holography. Digital hologram recording, numerical reconstruction, and related techniques, (Springer Berlin Heidelberg, 2010).

4. P. Ferraro, S. Grilli, D. Alfieri, S. De Nicola, A. Finizio, G. Pierattini, B. Javidi, G. Coppola, and V. Striano, “Extended focused image in microscopy by digital Holography,” Opt. Express 13(18), 6738–6749 (2005). [CrossRef] [PubMed]

5. J. Garcia-Sucerquia, W. B. Xu, S. K. Jericho, P. Klages, M. H. Jericho, and H. J. Kreuzer, “Digital in-line holographic microscopy,” Appl. Opt. 45(5), 836–850 (2006). [CrossRef] [PubMed]

6. D. K. Singh and P. K. Panigrahi, “Improved digital holographic reconstruction algorithm for depth error reduction and elimination of out-of-focus particles,” Opt. Express 18(3), 2426–2448 (2010). [CrossRef] [PubMed]

7. J. G. McNally, T. Karpova, J. Cooper, and J. A. Conchello, “Three-dimensional imaging by deconvolution microscopy,” Methods 19(3), 373–385 (1999). [CrossRef] [PubMed]

8. P. J. Shaw, “Comparison of wide-field/deconvolution and confocal microscopy for 3D imaging,” in Handbook of biological confocal microscopy, J. B. Pawley, ed., (Plenum Press, New York and London, 1995).

9. J. G. McNally, T. Karpova, J. Cooper, and J. A. Conchello, “Three-dimensional imaging by deconvolution microscopy,” Methods 19(3), 373–385 (1999). [CrossRef] [PubMed]

10. P. Sarder and A. Nehorai, “Deconvolution methods for 3-D fluorescence microscopy images,” IEEE Signal Process. Mag. 23(3), 32–45 (2006). [CrossRef]

11. W. Wallace, L. H. Schaefer, and J. R. Swedlow, “A workingperson’s guide to deconvolution in light microscopy,” Biotechniques 31(5), 1076–1097 (2001). [PubMed]

12. I. T. Young, “Image fidelity: Characterizing the imaging transfer function,” Methods Cell. Biol. 30, 2–47 (1989).

13. T. M. Kreis, “Frequency analysis of digital holography with reconstruction by convolution,” Opt. Eng. 41(8), 1829–1839 (2002). [CrossRef]

14. H. Meng and F. Hussain, “Holographic particle velocimetry: a 3D measurement technique for vortex interactions, coherent structures and turbulence,” Fluid Dyn. Res. 8(1-4), 33–52 (1991). [CrossRef]

15. S. Murata and N. Yasuda, “Potential of digital holography in particle measurement,” Opt. Laser Technol. 32(7-8), 567–574 (2000). [CrossRef]

16. R. Gold, “An iterative unfolding method for matrices,” Tech. Rep. ANL-6984 Argonne National Laboratory, Argonne, Illinois.

17. J. W. Goodman, “Introduction to Fourier optics,” (Roberts & Company Publishers, Engelwood, Colorado) 2005, Chap. 3.

18. S. A. Saq'an, A. S. Ayesh, A. M. Zihlif, E. Martuscelli, and G. Ragosta, “Physical properties of polystyrene/alum composites,” Polym. Test. 23(7), 739–745 (2004). [CrossRef]

Depth-resolved holographic reconstructions by three-dimensional deconvolution

Abstract

Corrections

1. Introduction

2. Image formation as convolution with the point-spread function

3. Formation and reconstruction of a hologram

4. Three-dimensional deconvolution of the reconstructed optical fields

4.1 Instant 3d-deconvolution of the reconstructed intensities by the Wiener filter

4.2 Iterative deconvolution of the reconstructed complex fields

5. Numerical simulation and reconstruction of holograms

6. Three-dimensional deconvolution of simulated holograms

6.1 Instant 3d-deconvolution of the reconstructed intensities by the Wiener filter

6.2 Iterative 3d-deconvolution of the reconstructed complex field

7. Experimental optical holograms

7.1 Instant 3d-deconvolution of the reconstructed intensities using the Wiener filter

7.2 Iterative 3d-deconvolution of the reconstructed complex fields

7.3 Three-dimensional deconvolution of the reconstructed optical fields with a simulated PSF

7.4 Iterative 3d-deconvolution of the reconstructed complex fields applied to extended objects

8. Conclusions

Acknowledgements

References and links

Cited By

Figures (12)

Equations (36)

Optics Express