Improved particle position accuracy from off-axis holograms using a Chebyshev model

Johan Öhman; Mikael Sjödahl

doi:10.1364/AO.57.00A157

1. INTRODUCTION

Digital holography is a suitable technique for tracking micro-particles since it records 3D information. Holography was originally proposed by Gabor in 1948 [1], and its digital version is usually attributed to Goodman [2]. Extracting information about particle fields is a natural application of digital holography because 3D information about position, size, and morphology is recorded in the hologram. Applications include studies of flow phenomena using holographic particle velocimetry (HPIV) [3], tracking microorganisms [4], and characterizing aerosols [5–7] and the dynamic behavior of combustion fuel particles [8,9] and air bubbles [10], to mention a few. There are two main configurations of holographic setups used for particle analysis. In in-line holography, the illumination wave is also used as a reference wave. Here the scattered light from particles interferes with the un-scattered light in the forward direction. This is a very simple and compact setup, but it has disadvantages in that only the forward direction can be used, and a twin image and zero-order terms disturb the reconstruction. The other main configuration is off-axis holography, which uses a separate reference wave with an off-axis angle. This setup is more complex, but with it, it is possible to remove the twin image and zero-order content from the image of interest. Any imaging direction can be used with this setup, not only the forward direction. In this investigation, the off-axis configuration is utilized because it allows for a general registration direction. In digital holography, the hologram is recorded on a CCD or CMOS detector and numerically reconstructed to obtain 3D information [11,12]. For low numerical apertures, the resolution in the reconstruction along the optical axis is poor compared to that of the lateral one [13]; therefore, particles appear elongated. To find the correct axial particle position, a metric of some kind is needed. Previously reported metrics are based on intensity [14–16], which are most common, or they are based on the whole complex field [17,18]. It is also possible to solve the inverse problem directly from the recorded hologram [19,20]. This is, however, more time consuming and, therefore, unsuitable for large particle sets.

We have, in a previous paper, presented a metric that is based on the wavefront curvature of the scattered light [21]. This metric utilizes the fact that scattered light from a particle converges toward its true position and thereafter diverges. The change from a converging to a diverging wave is found by tracking the wavefront curvature. One significant advantage with a phase-based metric is that it does not have an offset that is present in intensity-based metrics [22]. The wavefront curvature in our previous paper was estimated from phase gradients that were computed using finite differences. This was a simple way to estimate the curvature; unfortunately, in the presence of noise, the precision decreases rapidly due to the noise enhancement. In this paper, we propose a method that parameterizes the phase anomaly compared to a plane wave, without using finite differences. The method consists of fitting the phase anomaly to a three-dimensional polynomial model, and from the polynomial coefficients the particle position is determined. By using a model, we reduce the effect of noise. The rest of the paper is structured as follows: in Section 2, theory regarding particle scattering, digital hologram formation, and reconstruction are presented. The proposed polynomial model is also outlined. In Section 3, the setup is described along with the simulations and measurements performed. In Section 4, the results from both simulations and measurements are presented and discussed. Finally, the paper is summarized in Section 5, and the conclusions are presented.

Portions of this work were presented at the Digital Holography & 3-D Imaging conference in 2017, Paper M4A.3.

2. THEORY

A. Particle Holograms

Light incident on a particle generates an angularly dependent outwards propagating scattered field. Consider the geometry sketched in Fig. 1. The particle is located at the origin of a right handed Cartesian coordinate system and is illuminated by a plane wave propagating along the $y$ -axis. Detection is done along the vicinity of the $z$ -axis, which is denoted as the optical axis. We denote the Jones vectors of the incident light ${\bar{E}}_{i}$ and the scattered light ${\bar{E}}_{s}$ . For spherical micrometer-sized particles, the Lorenz-Mie theory is a suitable scattering model [23]. The theory is most commonly expressed with respect to the scattering plane; this is the plane spanned by the illumination and observation directions, denoted ${\bar{s}}_{i}$ and ${\bar{s}}_{s}$ , respectively, see Fig. 1. The incident and scattered fields are then expressed in polarization components parallel and perpendicular to the scattering plane. These components are denoted $E_{i ⊥}$ , $E_{i ∥}$ for the incident light and $E_{s ⊥}$ , $E_{s ∥}$ for the scattered light, respectively. The relationship between these components is as follows:

{\bar{E}}_{s} = (\begin{matrix} E_{s ⊥} (θ) \\ E_{s ∥} (θ) \end{matrix}) = \frac{e^{i k (r - z)}}{- i k r} (\begin{matrix} S_{1} (θ) & 0 \\ 0 & S_{2} (θ) \end{matrix}) (\begin{matrix} E_{i ⊥} \\ E_{i ∥} \end{matrix}),

where

S_{1} (θ)

and

S_{2} (θ)

are the angular dependent scattering coefficients,

θ

is the in-plane angle,

z

is the position along the optical axis, and

r = \sqrt{x^{2} + y^{2} + z^{2}}

is the distance from the particle. For the derivation of Eq. (1) and the full expression for the scattering coefficients, we refer the reader to [23].

Fig. 1. Scattering geometry with the defined scattering plane, spanned by vectors $s_{i}$ and $s_{s}$ . Polarization is expressed on component form, parallel, and orthogonal to this plane.

Download Full Size | PDF

In digital particle holography, the object light $\bar{E}$ is given by Eq. (1), and a reference wave is given by

\bar{R} = (\begin{matrix} R_{∥} \\ R_{⊥} \end{matrix}),

the interfere on the detector to form a digital hologram. Only components with mutual polarization interfere. In an off-axis holographic setup the object light, in this case the light scattered from particles, and the reference light are assumed to be separated in the angular spectrum [11]. In our setup, the light is recorded by the two-lens telecentric imaging system, depicted in Fig. 2. The magnification

M

between the object and image space is then given by

M = f_{2} / f_{1}

, where

f_{1}

and

f_{2}

are the focal lengths of lens

L_{1}

and lens

L_{2}

, respectively. We choose to denote the coordinate system in the object space by

(x, y, z)

and the coordinates in the image plane on the detector by

(x^{'}, y^{'}, z^{'})

. The imaging system then maps from the object plane to the image plane. The coordinates are related by the magnification

M

of the system as

x^{'} = - M x

,

y^{'} = - M y

, and

Δ z^{'} = M^{2} Δ z

, where

Δ z

and

Δ z^{'}

are the defocus distances in object and image space, respectively. The reference wave is an off-axis point source in the aperture plane that is reshaped by the lens

L_{2}

to become a plane wave. In these notations, the intensity on the detector for an arbitrary polarization becomes

I (x^{'}, y^{'}) = {| \bar{E} (x^{'}, y^{'}) \bar{R} (x^{'}, y^{'}) |}^{2} = {| E_{∥} (x^{'}, y^{'}) |}^{2} + {| E_{⊥} (x^{'}, y^{'}) |}^{2} + {| R_{∥} (x^{'}, y^{'}) |}^{2} + {| R_{⊥} (x^{'}, y^{'}) |}^{2} + E_{∥} {(x^{'}, y^{'})}^{*} R_{∥} (x^{'}, y^{'}) + E_{⊥} {(x^{'}, y^{'})}^{*} R_{⊥} (x^{'}, y^{'}) + E_{∥} (x^{'}, y^{'}) R_{∥}^{*} (x^{'}, y^{'}) + E_{⊥} (x^{'}, y^{'}) R_{⊥}^{*} (x^{'}, y^{'}),

where the first four terms are the intensity of each polarization component of the object wave and reference wave, respectively. The fifth and sixth terms are interference between the conjugate object wave and reference wave, for each polarization, respectively. Likewise, the seventh and eighth terms are the interference terms between the object wave and conjugate reference wave, for each polarization. By controlling the off-axis angle of the reference wave and the aperture size in the system, one can ensure no overlap between interference terms in the Fourier domain. Figure 3(a) shows the recorded intensity on the detector for a simulated hologram. The particle is defocused at distance

Δ z = 2000 μm

.

Fig. 2. Telecentric imaging system consisting of lenses $L_{1}$ and $L_{2}$ and aperture $A$ . Solid rays are the object light, and dashed rays are the reference light. Dashed black lines indicate the detector in image space and corresponding conjugate plane in object space.

Download Full Size | PDF

Fig. 3. Different properties in the reconstruction. (a) Intensity on the detector for a simulated hologram, (b) reconstructed phase along the center of the particle, (c) corresponding phase anomaly to (b). (d) Wavefront curvature at $Δ z = - 100 μm$ (black), $- 50 μm$ (blue), 0 μm (red), 50 μm (green), and 100 μm (magenta). In both (c) and (d), solid lines are polynomial fit, and markers the simulated data.

Download Full Size | PDF

B. Reconstruction

The hologram captured by the detector has to be reconstructed before further analysis is performed. It is possible to extract individual interference terms when using off-axis holography, since they are separated into different lobes in the Fourier domain. The first four terms in Eq. (3), that is, the intensity components, are located in a central lobe. The four interference terms each have a separate side lobe in a well configured system. The interference term of interest is extracted by applying a rectangular mask around the corresponding lobe, removing all content outside of it. If the reference wave $\bar{R}$ is a tilted plane wave, it will scale the object wave and add a phase gradient proportional to the off-axis angle. It is possible to remove the phase gradient by re-centering the lobe after the filtering process. By doing so, $R$ reduces to a scalar depending on the reference wave intensity. Because of this, we drop the $R$ from now on; hence, the complex amplitude $\tilde{E} (f_{x}^{'}, f_{y}^{'}; z_{D}^{'})$ is known in Fourier space for the component of interest, where $Δ z_{D}^{'}$ is the axial position of the detector in image space. It is now possible to compute the field at an axial position $z^{'} = z_{D}^{'} + Δ z^{'}$ from the detector with the use of the angular spectrum method [11], and the field becomes

E (x^{'}, y^{'}, z^{'}) = F^{- 1} [\tilde{E} (f_{x}^{'}, f_{y}^{'}; z_{D}^{'}) \exp (- j k s_{z}^{'} Δ z^{'})],

where

F^{- 1}

is the 2D inverse Fourier transform,

k = \frac{2 π}{λ}

is the wave number, and

s_{z}^{'} = \sqrt{1 - {(λ f_{x}^{'})}^{2} - {(λ f_{y}^{'})}^{2}}

is the directional cosine component in the z-direction. By reconstructing a stack of images separated by a distance

Δ z^{'}

, a discrete volume is produced. This reconstruction process creates a 3D volume from a 2D hologram. The volume contains information about both the amplitude and phase from which the 3D particle positions can be found using a suitable metric.

The basic principle of our metric is that the reconstructed wave from a particle converges toward its position and thereafter diverges. Converging and diverging waves have different signs of curvature, so the plane where the curvature changes sign is then the plane where the particle is located. The wavefront curvature can be obtained from the phase in the reconstructed volume. The phase in the volume oscillates quickly along the optical axis, much faster than the step length $Δ z^{'}$ in the reconstruction. This is problematic, since it then contains discontinues in the z-direction, making it unsuitable to fit it to a polynomial model. Figure 3(b) shows the reconstructed phase along the z-direction along the center of the particle. Aliasing effects are present, causing the oscillations to appear slower than they really are. By instead comparing the phase to a plane wave situated in the same plane, it becomes possible to compare the data between two planes. This phase anomaly is calculated as follows:

ϕ_{a} (x^{'}, y^{'}, z^{'}) = ∠ [E (x^{'}, y^{'}, z^{'})] - k z^{'},

where

k

is the wave number. We will from now on use the phase anomaly

ϕ_{a}

instead of the pure phase, since it is a smoother function more suitable for fitting. We now have two sets of information, the 3D intensity distribution

I (x^{'}, y^{'}, z^{'}) = {| E (x^{'}, y^{'}, z^{'}) |}^{2}

, and the 3D phase anomaly

ϕ_{a} (x^{'}, y^{'}, z^{'})

, from which the particle position can be estimated.

C. Chebyshev Polynomial Model

It is assumed that the phase anomaly around a single particle only depends on the response from that particle. Only a small volume around a potential particle is therefore used in the parameterization. This assumption holds in general if the particle distribution is sparse. Potential particles are found from a coarse intensity-based search where voxels with an intensity over a defined threshold are chosen.

For a single particle the phase anomaly, $ϕ_{a}$ , is fit to a model based on Chebyshev polynomials of the first kind [24]. Chebyshev polynomials are chosen as base functions for the parameterization for two reasons. The first is that they are orthonormal in a Cartesian coordinate system and, therefore, suitable for computations in a rectangular volume. Second, terms of different orders are associated with properties in the imaging process when using small to moderate numerical apertures. For instance, first-order terms are related to a spatial shift in direction, and second-order terms are related to defocus. The modeled phase anomaly $\hat{ϕ}$ in three dimensions is defined from one-dimensional Chebyshev polynomials. The polynomials are defined in the interval $[- 1, 1]$ , and the coordinates $[x^{'}, y^{'}, z^{'}]$ are rescaled to normalized coordinates $[ϵ, η, ζ]$ in this domain. The 3D polynomial is then as follows:

{\hat{ϕ}}_{a} (ε, η, ζ) = \sum_{i = 0}^{n} \sum_{j = 0}^{n} \sum_{k = 0}^{n} a_{i j k} T_{i} (ε) T_{j} (η) T_{k} (ζ),

where

T_{i}

,

T_{j}

, and

T_{k}

are Chebyshev polynomials of order

i

,

j

,

k

, respectively, and

a_{i j k}

are the associated coefficients for each base function triplet. Equation (6) can be rewritten in matrix form as

{\hat{ϕ}}_{a} (ε, η, ζ) = T a .

The objective is then to find the set of estimated coefficients ${\hat{a}}_{i j k}$ that best represents the measured phase anomaly $ϕ_{a}$ . We minimize the squared error between the measured data, $ϕ_{a}$ and the estimated data, ${\hat{ϕ}}_{a}$ , from the polynomial model. The coefficients are then the solution to the following minimization problem:

{\hat{a}}_{i j k} = \underset{a}{argmin} \sum W {(ϕ_{a} - {\hat{ϕ}}_{a})}^{2},

where

W

are weights that influence how much each data point contributes to the estimation. We chose the intensity of the scattered light as these weights. This ensures that only data points where there is light influence the estimation. The solution to this minimization problem is

\hat{a} = {(T^{T} W T)}^{- 1} T^{T} W ϕ_{a},

which is most efficiently calculated using Gauss elimination. To ensure that the estimated coefficients,

{\hat{a}}_{i j k}

, in Eq. (6), represent the phase anomaly correctly, a sufficiently high polynomial degree

n

is needed. To find this value, scattered light from a single particle is simulated using Eqs. (1) and (3), respectively. For the full description of the simulation implementation, we refer the reader to [21]. The simulated hologram is reconstructed according to Section 3. In the reconstructed volume, the parameterization is performed and the mean squared error (MSE) per voxel is evaluated for different values of

n

. The result is shown in Fig. 4. A good fit is obtained for

n = 4

, where the mean error is approximately

10^{- 3} rad / voxel

. Using higher orders does not significantly improve the fit. Therefore, a fit with

n = 4

is assumed to be sufficient and will be used further on in this paper. The phase anomaly corresponding to Fig. 3(b) is shown in Fig. 3(c). The markers are the simulated data, and the solid line is the Chebyshev polynomial fit. Figure 3(d) shows phase anomaly cross-sections in the y-direction for different z-positions. Here the change in wavefront curvature becomes apparent.

Fig. 4. MSE per voxel for different orders of $n$ .

Download Full Size | PDF

The three-dimensional data is now represented by the coefficient array $a_{i j k}$ . The relationship between the elements in the array and information about the particle is now investigated. A single particle is simulated, and coefficients $\hat{a}$ are estimated for different axial positions by translating the computational box along the optical axis. The dimension of the box used is $5 \times 5 \times 21$ voxels. All coefficients are investigated, and as an example, the results for the coefficients $a_{200}$ and $a_{020}$ are shown in Fig. 5. The coefficient $a_{200}$ has a linear behavior around the particle and becomes zero at the particle position. $a_{200}$ is therefore used as a metric to find the axial position of the particle. This coefficient corresponds to the quadratic term in the direction perpendicular to the scattering plane, so it describes the curvature of the wavefront in that direction. When $a_{200}$ becomes zero, the wavefront changes curvature. This behavior is consistent with the findings in [21]. The corresponding coefficient $a_{020}$ , which represents the quadratic term in the direction parallel to the scattering plane, is also shown in Fig. 5. It shows a similar behavior, but it is slightly shifted compared to $a_{200}$ . Therefore, only $a_{200}$ will be used as a metric.

Fig. 5. Evaluation of square terms $a_{200}$ and $a_{020}$ around particle location.

Download Full Size | PDF

To use the metric on a set of particles recorded in a single hologram, some practical considerations have to be made. First, we need to segment particles from each other, which is done by applying a threshold on the intensity in the volume. From the binary volume, connected voxels are clustered and labeled. The point with maximum intensity has a slight offset from the particle position, but it is still in the vicinity [22]. This point is therefore a good starting point when searching for the plane where $a_{200}$ becomes zero. A small 3D domain is formed around the maximum intensity point, and the Chebyshev model is evaluated in this domain. The position where $a_{200}$ becomes zero is then the position of the particle. This position is found by computing the zero-crossing using linear interpolation. It is important that the particles are sparse, so that only scattered light from an individual particle appear in the 3D domain. The whole procedure for obtaining the axial position of particles from a hologram can therefore be summarized in the following steps:

1. Reconstruct a three-dimensional volume using the angular spectrum method.
2. Apply a threshold on the intensity in the volume to obtain a binary volume.
3. Automatic segmentation of each particle from the binary volume.
4. Find the point of maximum intensity for each particle.
5. Evaluate the coefficient $a_{200}$ in a small domain around the point of maximum intensity for each particle.
6. The axial location where the coefficient $a_{200}$ becomes zero is stored as the position of each particle.

3. EXPERIMENTS

A. Setup

Consider the setup shown in Fig. 6. A 532 nm Nd:YAG continuous laser with a power of 85 mW is used for illumination. The light is polarized parallel to the scattering plane. This polarization is chosen to reduce Rayleigh scattering from inhomogeneous parts in the sample. The light is split into two parts using a 50/50 beam splitter (BS). The first part of the light is expanded by the lenses L1 and L2, with focal lengths $f = 150 mm$ and $f = - 20 mm$ , respectively. The expanded light illuminates the sample, which consists of a $23 mm \times 23 mm \times 23 mm$ silicon cube with spherical particles molded into it. The particles have an average diameter of 10 μm, and the refractive indices of the silicon and particles are 1.4 and 1.5, respectively. The concentration of particles is low so that the imaged scene is sparse, and only a few particles are present. The cube is placed so that only scattered light from the particles is imaged. There could be internal variations of the refractive index in the silicone cube that affect the experimental results compared to if the particles were imaged in air. These will then show up as an increase in random error. Scattered light from the particles is recorded using a telecentric imaging system, formed by lenses L3 and L4 at a 90° angle. Lenses L3 and L4 have focal lengths $f = 60 mm$ and $f = 80 mm$ , respectively. This corresponds to a system magnification of $M = 1.33$ . A $3 mm \times 3 mm$ aperture is used to band limit the image. It is placed in the joint focal plane between lenses L3 and L4, making the system telecentric. The other light part from the BS is used as the reference wave. It is focused into a polarization maintaining fiber using a mirror (M) and a fiber port (FP). The fiber end is placed in the aperture plane. By doing so, we ensure that the reference light becomes a point source in the Fourier domain. In the back focal plane of lens L4, a detector (Sony XCL 5005) is placed to record the hologram. It has $2448 \times 2048$ pixels and a pixel pitch of 3.45 μm, giving the detector a total size of $8.48 \times 7.1 mm$ . The images are stored in an 8-bit format. The fiber has a numerical aperture of 0.12, which is sufficient to cover the whole detector. The numerical aperture of this system is ${NA}_{0} = 0.025$ in the object plane and ${NA}_{1} = 0.0187$ in the image plane. Since the sample is a silicon cube, the effective numerical aperture is even smaller, since the silicon-air boundary will refract the scattered light. The effective numerical becomes

{NA}_{0}^{'} = \frac{{NA}_{0}}{n_{silicon}} \approx 0.0179 .

Fig. 6. Experimental setup used in the recordings. BS is a beam splitter, L1 is a $f = 150 mm$ lens, L2 is a $f = - 20 mm$ lens, L3 is a $f = 60 mm$ lens, L4 is a $f = 80 mm$ lens, A is the aperture, M is a mirror, and FP is a fiber port.

Download Full Size | PDF

B. Simulations and Measurements

To compare how the parametric Chebyshev method compares to the previously reported finite difference method, a large set of simulations was performed. Once again we use the simulation routine described in [21], that is, an implementation of Eqs. (1) and (3), respectively. The simulation properties are set as described in Section 3.A. Some assumptions are made: the illumination is a perfect plane wave, and particles are identical perfect spheres with a diameter of 10 μm. Then, 10 holograms each containing 20 particles located at random lateral and axial positions are simulated. Since the axial positions are known, it is easy to compute the error for each particle. The standard deviation, false positive ratio, and false negative ratio are computed. This is repeated for different noise levels ranging from 0 to 4 percent of the maximum intensity.

For experimental measurements, the backside of the silicon cube is placed at approximately $L = 60 mm$ in front of the imaging system, so that the focus plane is inside the sample cube. Unlike in simulations, the true particle position is not known, and there are several holograms of stationary particles recorded in the series. Particle positions are evaluated with the method outlined in Section 4 from the holograms. A set of positions for five different particles are then obtained. For each particle, the axial mean position is estimated and subtracted from the recorded positions. The associated errors are therefore computed.

4. RESULTS AND DISCUSSION

The resulting standard deviations from the simulations are shown in Fig. 7 as a function of the noise ratio. Data for both the Chebyshev model and the finite difference method are presented. It is clear that by using the Chebyshev model, the standard deviation drastically decreases for larger noise ratios. At the highest simulated noise ratio of 4%, the standard deviation is reduced from 39 μm to 10 μm, an improvement of about 4 times. This can be explained by the fact that the finite difference is influenced by the noise in two voxels, while the Chebyshev polynomial method uses far more data points. Since the noise is assumed to be zero mean and independent, it will have a smaller effect with the use of more data points. For no or low noise ratios, the proposed polynomial model performed worse than the finite difference method. This is likely due to the extension from 2D to 3D computational domain. In a volume, it is more likely that light from other particles influences the estimation. The false positive ratio is around 1%, so it is very unlikely that noise is mistaken for a particle. This is due to the fact that noise rarely originates from a geometrical point, and hence, it does not have the behavior of an outwards propagating wavefront. Therefore, it does not have the same change in wavefront curvature as a particle.

Fig. 7. Simulated standard deviations of axial position error using both Chebyshev model (blue) and finite difference method (red).

Download Full Size | PDF

The false negative ratio increased to around 5% when using the Chebyshev model, compared to 3% when using finite difference. A closer study of the undetected particles shows two main reasons for the increase. The first reason is inaccuracy in the course estimation. The model searches for a zero-crossing of $a_{200}$ close to the point obtained from the course estimation. So if that estimation is very inaccurate, the model never finds the zero-crossing. The second reason is that particles located near the edges of the reconstructed volume do not use the same number of data points, since parts of the computational domain are shifted outside the reconstruction. This can cause the model to behave inaccurately.

The performance is also compared to the total sum of gradients method, which is commonly used to determine particles axial positions [25,26]. Applying the method to the simulated data set shows that the sum of the gradients method results in a small offset from the true axial position, around 3 μm, which is still inside the particle. The standard deviation for the different noise levels showed a slight increase ranging from 7–12 μm compared to 5–10 μm for the polynomial fit. The two methods, hence, show comparable noise susceptibility.

The results from the experimental measurements are presented in Fig. 8. This figure shows the error histogram after outliers have been removed. The error seems to have a Gaussian profile, where the standard deviation is $σ = 13 μm$ , which is an improvement compared to the finite difference method, where the standard deviation was 18 μm [21]. Even though the standard deviation decreased, it is still higher than for the simulations. From Fig. 7, the expected standard deviation would be around 7 μm. This can depend on several factors. In the simulations, we assume only single scattering and ignore any effects from the finite size of the silicon cube. For a single particle, scattered light from other particles can be viewed as noise. It is possible that these effects influence the results. Also, aspects like mechanical and thermal stability in the experimental setup influence the results. Still, the standard deviation is in the same order of magnitude as the diameter of the particle. It is, therefore, likely that the measured position is inside or in the close vicinity of the particle.

Fig. 8. Histogram of experimental error. Error is approximately zero mean with a standard deviation of 12.8 μm.

Download Full Size | PDF

Finally, a brief discussion of the parametrization process is in order. Consider the errors presented in Fig. 4. For $n = 4$ , the error is negligible. However, since the coefficient used as a metric is a second-order term, higher terms are redundant due to the orthonormal properties of the base functions. This means adding higher order terms should not change lower order coefficients. The biggest error reduction is from $n = 1$ to $n = 2$ . This is not a surprise, since a spherical wave over a small aperture can be approximated by a second-order term. This is also why a second-order term $a_{200}$ is used as a metric. From Fig. 5, it is clear that only $a_{200}$ can be used as a metric, since $a_{020}$ is shifted along the optical axis. These coefficients represent quadratic terms parallel and perpendicular to the scattering plane, respectively. The wavefront parallel to the scattering plane is modulated by Lorenz-Mie scattering coefficients, while the perpendicular direction is not. The modulation distorts the wavefront in the parallel direction, causing $a_{020}$ to become zero at a shifted position, while $a_{200}$ is unaffected. This effect has previously been reported by Pu and Meng [27]. They related the modulation effects to aberrations for different observation directions. These effects are important to consider when interpreting results from a certain observation direction. A future extension of the proposed polynomial method to include non-spherical particles could be possible. Non-spherical particles do, in general, depend on the perpendicular direction as well [23], so the modulation effects need to be handled. This is something that could be addressed in future work. Modulation effects can also be introduced by the reference wave. If the reference is a curved wave instead of a plane wave, the curvature will modulate the phase response. It is therefore important to ensure the quality of the reference wave.

5. CONCLUSIONS

In this paper, a parameterized model based on Chebyshev polynomials is introduced to estimate the wavefront curvature. It is shown that the phase anomaly around a scattering particle can be modeled by a three-dimensional Chebyshev polynomial of degree $n = 4$ . To find this degree, the mean error per voxel is evaluated. The polynomial coefficients hold information about the particle position, and by using simulations, the behavior of all obtained coefficients was studied. It is found that the coefficient $a_{200}$ becomes zero in the axial plane where the particle is located. This coefficient corresponds to the quadratic term perpendicular to the scattering plane. By evaluating this coefficient at different positions and searching for the axial plane where it becomes zero, it is possible to find the true particle position. From simulations, this method has standard deviations between 6 and 10 μm. This should be compared to the standard deviation using the finite difference method previously reported, which were 3–39 μm for the same conditions. The Chebyshev polynomial model performs better in higher noise ratios, since the noise in individual voxels does not affect the result to the same degree. For no or very low noise, the worse performance is due to the extension of the computational domain from 2D to 3D. In a volume, light from other particles is more easily introduced in the computation. Experimental results also show improvement when using this method, from 18 μm to 13 μm. This improvement is not as large as the improvement expected from simulations; however, ideal conditions in simulations are rarely obtainable.

Funding

Vetenskapsrådet (VR) (621-2014-4906).

REFERENCES

1. D. Gabor, “A new microscopic principle,” Nature 161, 777–778 (1948). [CrossRef]

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

3. H. Meng, G. Pan, Y. Pu, and S. H. Woodward, “Holographic particle image velocimetry: from film to digital recording,” Meas. Sci. Technol. 15, 673–685 (2004). [CrossRef]

4. J. Sheng, E. Malkiel, and J. Katz, “Digital holographic microscope for measuring three-dimensional particle distributions and motions,” Appl. Opt. 45, 3893–3901 (2006). [CrossRef]

5. M. J. Berg and S. Holler, “Simultaneous holographic imaging and light-scattering pattern measurement of individual microparticles,” Opt. Lett. 41, 3363–3366 (2016). [CrossRef]

6. M. J. Berg, N. R. Subedi, and P. A. Anderson, “Measuring extinction with digital holography: nonspherical particles and experimental validation,” Opt. Lett. 42, 1011–1014 (2017). [CrossRef]

7. O. Kemppinen, Y. Heinson, and M. Berg, “Quasi-three-dimensional particle imaging with digital holography,” Appl. Opt. 56, F53–F60 (2017). [CrossRef]

8. Y. Wu, M. Brunel, R. Li, L. Lan, W. Ao, J. Chen, X. Wu, and G. Gréhan, “Simultaneous amplitude and phase contrast imaging of burning fuel particle and flame with digital inline holography: model and verification,” J. Quant. Spectrosc. Radiat. Transfer 199, 26–35 (2017). [CrossRef]

9. Y. Wu, L. Yao, X. Wu, J. Chen, G. Gréhan, and K. Cen, “3D imaging of individual burning char and volatile plume in a pulverized coal flame with digital inline holography,” Fuel 206, 429–436 (2017). [CrossRef]

10. S. I. Satake, Y. Yonemoto, T. Kikuchi, and T. Kunugi, “Detection of microbubble position by a digital hologram,” Appl. Opt. 50, 5999–6005 (2011). [CrossRef]

11. J. W. Goodman, Introduction to Fourier Optics, 2nd ed. (McGraw-Hill, 1996).

12. N. Verrier and M. Atlan, “Off-axis digital hologram reconstruction: some practical considerations,” Appl. Opt. 50, H136–H146 (2011). [CrossRef]

13. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light, 7th ed. (Cambridge University, 1999).

14. Y. S. Bae, J. I. Song, and D. Y. Kim, “Volumetric reconstruction of Brownian motion of a micrometer-size bead in water,” Opt. Commun. 309, 291–297 (2013). [CrossRef]

15. J. K. Abrantes, M. Stanislas, S. Coudert, and L. F. A. Azevedo, “Digital microscopic holography for micrometer particles in air,” Appl. Opt. 52, A397–A409 (2013). [CrossRef]

16. L. Wilson and R. Zhang, “3D localization of weak scatterers in digital holographic microscopy using Rayleigh-Sommerfeld back-propagation,” Opt. Express 20, 16735–16744 (2012). [CrossRef]

17. G. Pan and H. Meng, “Digital holography of particle fields: reconstruction by use of complex amplitude,” Appl. Opt. 42, 827–833 (2003). [CrossRef]

18. J. de Jong and H. Meng, “Digital holographic particle validation via complex wave,” Appl. Opt. 46, 7652–7661 (2007). [CrossRef]

19. F. Soulez, L. Denis, C. Fournier, É. Thiébaut, and C. Goepfert, “Inverse-problem approach for particle digital holography: accurate location based on local optimization,” J. Opt. Soc. Am. A 24, 1164–1171 (2007). [CrossRef]

20. J. Fung, K. E. Martin, R. W. Perry, D. M. Kaz, R. McGorty, and V. N. Manoharan, “Measuring translational, rotational, and vibrational dynamics in colloids with digital holographic microscopy,” Opt. Express 19, 8051–8065 (2011). [CrossRef]

21. J. Öhman and M. Sjödahl, “Off-axis digital holographic particle positioning based on polarization-sensitive wavefront curvature estimation,” Appl. Opt. 55, 7503–7510 (2016). [CrossRef]

22. F. C. Cheong, B. J. Krishnatreya, and D. G. Grier, “Strategies for three-dimensional particle tracking with holographic video microscopy,” Opt. Express 18, 13563–13573 (2010). [CrossRef]

23. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley, 1998).

24. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (Dover, 1972), Vol. 9.

25. Y.-S. Choi, K.-W. Seo, M.-H. Sohn, and S.-J. Lee, “Advances in digital holographic micro-PTV for analyzing microscale flows,” Opt. Lasers Eng. 50, 39–45 (2012). [CrossRef]

26. P. Memmolo, L. Miccio, M. Paturzo, G. D. Caprio, G. Coppola, P. A. Netti, and P. Ferraro, “Recent advances in holographic 3D particle tracking,” Adv. Opt. Photon. 7, 713–755 (2015). [CrossRef]

27. Y. Pu and H. Meng, “Intrinsic aberrations due to Mie scattering in particle holography,” J. Opt. Soc. Am. A 20, 1920–1932 (2003). [CrossRef]

Improved particle position accuracy from off-axis holograms using a Chebyshev model

Abstract

1. INTRODUCTION

2. THEORY

A. Particle Holograms

B. Reconstruction

C. Chebyshev Polynomial Model

3. EXPERIMENTS

A. Setup

B. Simulations and Measurements

4. RESULTS AND DISCUSSION

5. CONCLUSIONS

Funding

REFERENCES

Cited By

Figures (8)

Equations (10)

Applied Optics