1. INTRODUCTION

In traditional digital particle image velocimetry, particle positions are only obtained within a 2D plane [1]. By utilizing holographic recordings, particle positioning can be done in a volume which opens up possibilities for investigatigating true 3D phenomenas [2]. This is possible since holographic imaging records the full characteristics of the light, i.e., both the amplitude and phase [3]. The most commonly used holographic setup for particle imaging is the in-line setup, where the illumination and reference beams are coaxial. This enables a very simple setup where no adjustments to the reference beam are necessary. The drawback with this setup is that not just the interference term between the object and reference light is present on the sensor. The directly transmitted intensity is also recorded. The twin image problem is also present in the hologram. To avoid these problems, off-axis holography can be used. In off-axis holography a portion of the illumination light is used as a separate reference beam that is incident on the sensor at an angle. By using this method, the interference term between the reference and object light can be filtered out and the other terms are removed. By recording the hologram on a CCD detector, it is possible to numerically reconstruct the light field at a distance away from the focus plane and hence create a reconstructed volume around the imaged scene [4]. However, the axial resolution is significantly worse than the lateral one at least eight times and scales with the numerical aperture [5]. This causes spherical particles to be cigar-shaped. The actual axial location of the particle is therefore not immediately obtained from the reconstruction. The problem of poor axial resolution in holographic reconstructions may seem a bit counterintuitive, as one application of traditional interferometry is measurement of height deviation over smooth surfaces with a precision of fractions of a wavelength. These traditional surface measurements, however, detect relative variation in a distance over a surface, whereas the absolute position of the surface remains unknown.

Locating particles in 3D is a well studied subject and several different methods have been proposed. The most common metric used to determine axial position is based on the intensity of the reconstructed light. To separate particles in the recorded data, a threshold process is applied to the intensity of the reconstructed volume. Only voxels with an intensity higher than the threshold will be considered to contain useful data. From this process, a binary volume is obtained. Using the binary volume, particles can be segmented and the data for each particle is labeled. The simplest metrics only utilize the intensity in the reconstructed volume. To improve resolution, the intensity of the lateral centroid or the overall intensity can be evaluated for each particle [6,7]. One can also utilize the intensity gradient to determine where the particle is at focus [7,8]. If it is known that the particles of interest are spherical, then the scattered light can be modeled with the Lorenz–Mie theory and a parameter fit of position, refractive index, and size can be done [9]. However, a general particle cannot be modeled by the Lorenz–Mie theory and such a parameter fit will only be approximate. A problem that arises when using an intensity metric is that the estimated axial position does not coincide with the true axial position. Cheong et al. studied this problem and compared the centroid of the reconstructed intensity to a fit of the exact Lorenz–Mie model for spherical particles [10]. The centroid position coincides with the position from the Lorenz–Mie fit for the lateral coordinates. For the axial position, the centroid had an offset and therefore would not yield the true position. It was found that the offset depends on the size of the particle and increases with increasing particle radius. This indicates that the offset could be used for particle sizing. By calibration, this offset can be removed in order to obtain the true position. One can also overcome this problem by making tomographic recordings using two perpendicular cameras. By combining the two reconstructed volumes, the axial resolution can be reduced to the lateral one [11].

The aim of this paper is to directly obtain the true axial position using the reconstructed phase of the scattered light. There have been few studies on particle positioning that is based on phase measurement. The behavior of the phase around the focus plane of a particle was studied in detail in previous studies for in-line holography [12] and complex metrics have been used to determine particle locations in in-line holography for both forward and side scattering setups [13,14]. The idea outlined in this paper is that the wavefront curvature holds information about the axial position of particles and that the true position can be estimated from the reconstructed phase.

This paper is structured in the following way; in Section 2 a simulation model for particle scattering is presented. In Section 3 the proposed metric for estimating particle locations is outlined. Simulations are done to show the behavior of the metric; results are presented and discussed. In Section 4 the experimental setup used to record holograms is outlined, and in Section 5 results from experimental measurements are presented and discussed. The localization metric is also applied to track an axial translation of a particle sample. Finally, in Section 6, conclusions are made based on the results from this paper.

2. THEORY

A. Simulation of Mie Scattering

Using numerical simulations, the behavior of the scattered light can be studied in detail. Scattering from spherical micrometer-sized particles can be modeled by Mie scattering. We will briefly cover this topic; for a more complete description we refer the reader to [15]. For a plane wave incident on a particle, the scattered light will consist of an outgoing wave modulated by Mie coefficients. The majority of the scattered light is concentrated to a narrow cone in the forward scattering direction. However, when imaging scattered light in the forward direction, much stronger illumination light is also present on the sensor. The consequence of this is that the shadow of the particle is recorded instead of the scattered light from the particle. Furthermore, the narrow light cone limits the spatial frequencies available for imaging. This means that the spatial frequencies are limited regardless of the properties of the imaging system. To overcome both these problems, imaging can be performed at a 90° angle. At this angle, the illumination light is not present on the sensor and it is not concentrated to a narrow cone. This light is, however, much weaker than the forward scattered light, which requires longer exposure times or higher illumination intensity.

According to the Lorenz–Mie theory, the scattered light is related to the incident light by the coefficients $S_1(\theta )$ and $S_2(\theta )$. These angularly dependent scattering coefficients determine how the scattered light from a spherical particle behaves in a certain direction. The parameter $S_1(\theta )$ relates how the polarization component orthogonal to the scattering plane behaves, and the parameter $S_2(\theta )$ relates how the component parallel to the scattering plane behaves. The scattering plane is defined as the plane spanned by the vector $s_i$ in the illumination direction and the vector $s_s$ in the scattering direction; see Fig. 1. The angle $\theta$ is the angle between these two vectors. The parallel and orthogonal components of the incident light are also indicated. $S_1(\theta )$ and $S_2(\theta )$ are, in general, complex variables so the scattering coefficients will affect both the amplitude and phase. The relation between the incident and scattered light is given by

 

Fig. 1. Definition of scattering plane with the illumination direction $s_i$ and scattering direction $s_s$.

Download Full Size | PDF

Assuming the incident light $E_i$ to be a plane wave linearly polarized parallel to the scattering plane, $E_i=E_{i\Vert }$, the scattered light at an angle $\theta$ is then given by

The set of scattering angles is centered around the optical axis. For the 90° scattering setup. the optical axis is along the $z$ axis; see Figs. 1 and 2(a). In general, the particle is not located at the focus plane, and therefore a defocus term is also present. Only scalar waves are used in this paper, so we only use one component of the scattered light and denote it by $U$ from now on. This simplification is reasonable, since the detection can be approximated to be done in the scattering plane. Considering the imaging geometry sketched in Fig. 2(a), the field at the joint focal plane $U_a$ can be expressed as follows

 

Fig. 2. Coordinate systems, where $U_i$ is the field at the focus plane, $U_a$ is the field in the joint focal plane, and $U_o$ is the field at the detector. $L_1$ and $L_2$ are imaging lenses with focal lengths $f_1$ and $f_2$, respectively. $\mathrm{\Delta }z$ is the defocus distance of the particle. (a) Shows the coordinates in the object- and image-domains and (b) shows the propagation of a image plane by a distance $\mathrm{\Delta }{z}^{\prime }$.

Download Full Size | PDF

For a telecentric imaging system, the input field at the front focus plane is related to the output field at the back focus plane as follows:

The spatial frequencies are also rescaled with the magnification. In addition, $U_o$ is low-pass filtered by a hard aperture in the joint focal plane $U_a$ between the lenses. Combining Eqs. (3), (4), and (6) describes the field at the detector given the position and particle size in object space. By only considering single scattering terms, the contribution of several particles can be superimposed in the output field.

B. Off-axis Holography and Reconstruction

In digital holography the object light $U$ and a reference wave $R$ interfere on the detector to form a hologram. In an off-axis holographic setup, the object light, in this case the scattered light from the particles, and the reference light are separated in the angular spectrum [4]. We choose to denote the coordinate system in the object plane by $(x,y,z)$ and the coordinates in the image plane on the detector by $(x^{\prime },y^{\prime },z^{\prime })$; see Fig. 2(a). The coordinates are related by the magnification $M$ of the system $x^{\prime }=-Mx$ and $y^{\prime }=-My$. In these notations the intensity on the detector becomes

The first two terms are the intensity of the object wave and reference wave, respectively, and the last two terms are interference terms between the object wave and reference wave. There is no overlap between the terms in the Fourier domain. The complex amplitude of the object wave modulated by the reference wave is obtained by filtering in the Fourier domain. A rectangular filter is applied to extract one of the interference terms. In a well-aligned imaging system, the reference is a plane wave that only scales the object light. One can then, without further loss of generality, write

Therefore, the recorded interference terms contain a full description of the complex amplitude of the detected object light. By generating the field at the detector using Eq. (6) and making it interfere with a plane wave, the simulated hologram is produced. Thus the detected interference term $U(x^{\prime },y^{\prime })R^{*}(x^{\prime },y^{\prime })$ is formed.

Reconstruction of the target volume follows the same path as the simulation. By knowing the complex amplitude of the light field, one can propagate it to a different plane along the optical axis. We utilized the angular spectrum method [4]. Before further calculations the tilt phase contribution from the off-axis reference light is compensated for by centering the spatial frequency content. If the field at the detector is denoted $\mathit{U}(x^{\prime },y^{\prime },0)$, then the field at a distance $\mathrm{\Delta }{z}^{\prime }$ away from the detector plane is given by

3. PARTICLE ESTIMATION

Because of the low numerical aperture, the axial resolution in the volume is poor so that spherical particles become elongated and cigar-shaped in the reconstruction. Therefore, a metric is required to find the true position of the particle. Since the scattered light is an outgoing wave originating from the particle, the wavefront curvature changes its direction of propagation at the particle’s location. In the wavefront reconstruction, the scattered light converges toward the particle and thereafter diverges when passing it. This change in curvature is the core phenomena in the estimation. To more accurately determine a particle’s location, we suggest a metric that is based on the phase gradients of the scattered light rather than the phase itself. The idea is that the light scattered from a particle will have a plane wavefront at the particle location. Away from the particle, the wavefront has a curvature that is directed toward the particle. So when going through the focus, the curvature changes sign. The proposed metric can be described in the following way. The wavefront curvature needs to be estimated. This is done from the reconstructed phase in the volume. The imaging system will image a limited part of the wavefront, and the phase around a particle can then be approximated with a second-order polynomial. Still, the phase can be wrapped, and a direct second-order fit fails in that scenario. To avoid this problem, the in plane phase gradients are used instead of the phase itself. This will also allow for a lower-order fit later on that is less complex and has better convergence. The phase gradients are defined as follows:

A. Simulated Holograms

To investigate the assumption that the phase gradients in fact become zero at the true axial position, the behavior of a single particle is simulated. The simulated system is set to resemble the experimental setup. The particle is set to have a diameter of 10 μm and a refractive index of 1.5, and the particle is non-absorbing. The surrounding medium is set to be silicon with a refractive index of 1.4. The illumination light has a wavelength of 532 nm and is linearly polarized parallel to the scattering plane. The two imaging lenses are set to have focal lengths of 60 mm and 80 mm, respectively; that gives a lateral magnification of $M=1.33$. A square aperture of $3\times 3\mathrm{mm}$ is used to band-limit the light. This corresponds to a numerical aperture in the object plane of ${\mathrm{NA}}_0=0.025$ and ${\mathrm{NA}}_1=0.0187$ in the image plane. The pixel size of the detector is $d=3.45\mathrm{\mu m}$. The particle is located at $(x,y,z)=(0,0,2500)\mathrm{\mu m}$, i.e., in the center of the object plane with a defocus of 2500 μm toward the imaging system. The field is discretized in both lateral coordinates corresponding to the detector size and intensity values ranging from 0 to 255. From the acquired hologram, a volume is reconstructed using Eq. (9). To evaluate the phase gradients from this volume, segmentation of voxels corresponding to the particle is done. To do this, a threshold is set at 10% of the maximum intensity to form a binary image of the volume. Voxels are clustered to form the regions of interest for the particle. The parameters $B_2$ and $A_1$ are calculated using the method outlined above. In Fig. 3 parameter values are presented along with the centroid intensity of the particle. Here one can clearly see the shift in curvature that we are interested in. For a perfectly spherical wave emitted from a point source, the phase gradients are spherically symmetrical, so the phase gradients both parallel and orthogonal become zero in the focus plane. For particles of finite size, Mie effects are present in the scattered light. The Mie coefficients modulate the behavior of the light depending on the angular direction $\theta$ in the scattering plane. This angular dependence varies with pixel position parallel to the scattering plane and not with the orthogonal position. This modulation with respect to parallel position will affect the phase gradients parallel to the scattering plane, while not affecting the gradients orthogonal to it. The parameter $A_1$ is therefore distorted and its zero-crossing is shifted away from the true particle position. In contrast, the parameter $B_2$ is well behaved and linear. Away from the particle this parameter also becomes distorted. This is because fewer points are available for the parameter fit. In Figs. 3(a) and 3(b), both the parallel and orthogonal components are shown. Figure 3(b) is just a zoomed-in version of Fig. 3(a). The parameters show exponential behavior as a function of $z$; they behave approximately linear around the axial location of the particle. The linear slope around the zero-crossing will increase with increasing numerical aperture ${\mathrm{NA}}_0$, since a larger part of the wavefront is used in the detection. So if all spatial frequencies are available in the detection, a discontinuity will be present at the particle location. Since the detection in Fig. 3 is band-limited, the slope is finite. From this simulation we can conclude that the parameter $B_2$ becomes zero in the vicinity of the true axial position, while the zero-crossing for the parallel parameter $A_1$ is shifted. Therefore, only phase gradients orthogonal to the scattering plane can be used to locate particles accurately. In our coordinate system, the parallel direction is in the $x$ direction and the orthogonal direction is in the $y$ direction. All this is just an effect of the scattering parameters only being sensitive to the in-plane angle.

 

Fig. 3. Evaluation from simulated data; (a) shows the phase-gradient tilt for parallel direction (dashed) and orthogonal direction (solid), (b) shows a zoomed-in version of (a), and (c) shows the centroid intensity.

Download Full Size | PDF

Further simulations of 40 holograms with 40 particles each are performed. The lateral positions of the particles are chosen at random, while the axial is kept at $z=2500\mathrm{\mu m}$. Simulations are repeated for different noise levels ranging from 0% to 4% of the maximum bit depth in the detector, with an increment of 0.5%. In our experimental measurements, the exact particle positions are unknown and a translation is estimated instead. Therefore, a translation is simulated as well. One hologram is acquired before a 100 μm translation and a hologram is recorded afterward. From these simulations, the parameter $B_2$ is evaluated for each particle and the axial position of the zero-crossing is saved as the particle’s axial location. Also, a fit of the Gaussian profile for the centroid intensity is performed to obtain the location of the peak. The results from the simulation with a noise level of 1% are presented as histograms in Fig. 4. From these histograms it becomes more clear that the intensity peak is shifted approximately 20 μm, while the positions acquired from the phase-gradient method are centered around the true axial position. By fitting a Gaussian distribution to the data, the mean and standard deviations are calculated. For the centroid intensity method, the mean is $-2519.6\mathrm{\mu m}$ with a standard deviation of 4.5 μm, while for the phase-gradient method, the mean is $-2499.7\mathrm{\mu m}$ with a standard deviation of 5.1 μm. From these simulations we can conclude that the phase-gradient-based method produces the correct axial position of the particle and that both methods have comparable standard deviation. To evaluate the results from an axial translation particle, positions before and after the translation are paired. The resulting mean translation is 100 μm, as expected, and the standard deviation depends on the noise level. The resulting standard deviations are presented in Table 1. The standard deviation appears to scale approximately linear with the noise ratio. One additional feature with this phase-gradient metric is that it rejects noise that shows up as particles when using only intensity. Noise does not have the phase behavior of a particle since it is not passing a true geometrical focus. A similar validation process was utilized by de Jong and Meng [14], where they posed a criteria on the real part of the complex field. Although our metrics differ from the requirement on potential particles described in [14], they are the same as [14] in the sense that they show typical particle characteristics. By enforcing this requirement, one can use a rather low threshold since non-particle voxel clusters will be removed by the use of the phase-gradient metric. This method is therefore insensitive to false particles generated from noise. False negative and false positive particle detections are quantified and are between 1.5% and 3% for different noise levels. Compared with quotas previously reported, the false negative values are about the same, while the false positive quota is lower [14].

 

Fig. 4. Histograms of simulated data with Gaussian distribution fits for both phase-gradient method (a) and centroid intensity method (b). The solid lines are the fitted Gaussian distributions.

Download Full Size | PDF

Table 1. Simulated Standard Deviation for Different Peak-to-Peak Noise Ratios

View Table

4. EXPERIMENTAL MEASUREMENTS

The experimental setup consists of the off-axis holographic telecentric imaging system, depicted in Fig. 5. For illumination, a 532 nm Nd:YAG continuous laser with a power of 85 mW is used. The polarization is adjusted with a half-wave plate (not shown in the figure) so that the light is linearly polarized parallel to the scattering plane. Using a beam splitter (BS), part of the light is focused into a polarization-maintaining fiber using a mirror ($M$) and a fiber coupling (FC). The rest of the light is expanded by the biconcave lens L2 of focal length $f=-20\mathrm{mm}$ and the biconvex lens L1 of focal length $f=150\mathrm{mm}$, so that a larger volume of the sample is illuminated. The sample consists of a $23\mathrm{mm}\times 23\mathrm{mm}\times 23\mathrm{mm}$ silicon cube that has spherical particles molded into it. Refractive indices of the silicon and particles are 1.4 and 1.5, respectively, and the particles have an average diameter of 10 μm. The concentration of particles is low so that only a about 20–50 particles are present in the imaged scene. Imaging of the scattered light is performed at an angle of 90° using a telecentric imaging system. Lens L3 is a $f=60\mathrm{mm}$ biconvex lens and lens L4 is a $f=80\mathrm{mm}$ biconvex lens. This corresponds to a system magnification of $M=1.33$. In the common focal plane, a square aperture with the size of $3\mathrm{mm}\times 3\mathrm{mm}$ is placed to control the spatial frequency content in the image. By placing the aperture here, the entrance- and exit-pupils of the system are located at infinity, hence making the system telecentric, so that the magnifications are constant with respect to axial distance. To generate a hologram, the other end of the polarization-maintaining fiber, containing the reference light, is connected to the aperture plate. 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 CCD detector (Sony XCL 5005) is placed to record the hologram. The detector has $2448\times 2048\text{pixels}$ with a size of 3.45 μm, giving it an effective size of $8.48\mathrm{mm}\times 7.1\mathrm{mm}$. The fiber has a numerical aperture of 0.12, which is sufficient to cover the whole detector. The numerical aperture of this system is ${\mathrm{NA}}_0=0.025$ in the object plane and ${\mathrm{NA}}_1=0.0187$ in the image plane, respectively. 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 aperture becomes

 

Fig. 5. Experimental setup used in the recordings. BS is a beam splitter, L1 is a $f=150\mathrm{mm}$ lens, L2 is a $f=-20\mathrm{mm}$ lens, L3 is a $f=60\mathrm{mm}$ lens, L4 is a $f=80\mathrm{mm}$ lens, $A$ is the aperture, $M$ is a mirror, and FC is a fiber collimator.

Download Full Size | PDF

The backside of the silicon cube is placed at approximately $L=60\mathrm{mm}$ in front of the imaging system so that the focus plane is inside the sample cube and a hologram is recorded of the side-scattered light. It is then translated 100 μm along the optical axis away from the imaging system and another hologram is recorded.

5. EXPERIMENTAL RESULTS, AND DISCUSSION

From the holograms recorded, reconstructed volumes and phase gradients are estimated using the method outlined in Section 3. It is of interest if the parameter $B_2$ will coincide with the simulated behavior. The coefficient $B_2$ is evaluated for 17 particles and a new variable $z^{*}$ is introduced to compare individual particles regardless of axial position. It is defined as

 

Fig. 6. Evaluation of parameter $B_2$ for 17 particles from an experimental recording. The $z$ axis is shifted so that all particles have their zero-crossing at $z^{*}=0$.

Download Full Size | PDF

 

Fig. 7. 3D scatter plot of particle positions before and after the axial translation. Links between paired particles are also indicated.

Download Full Size | PDF

The estimated mean translation between recordings is 105.0 μm with a standard deviation of 25.3 μm. The standard deviation of the individual positioning is hence $25.3\mathrm{\mu m}/\sqrt{2}=17.9\mathrm{\mu m}$, assuming that the translation is the difference between two independent distributions. This indicates that the mean translation is in-line with the translation performed by the micrometer screw; the standard deviation is, however, quite significant. Comparing the estimated translations with the simulated translations, it is found that both are unbiased, so the mean translation is close to 100 μm. The standard deviation from the experimental data is in the same order of magnitude as the simulated values in Table 1. The experimental standard deviation corresponds to a simulated noise level of approximately 2.5%–3%. From this we can conclude that the noise introduces an uncertainty in the phase-gradient estimations that causes an uncertainty in the estimated position. To improve the setup, a larger numerical aperture can be used. Then, a larger portion of the wavefront is imaged, and hence the curvature would be clearer. A larger magnification will also make the positioning more accurate, since the image and object domains are related by the magnification.

One important aspect of the phase-gradient method is that it rejects potential particles that do not have the particle signature. This will reject potential particles that originate from noise, and hence the number of false particle detections is reduced. Actual particles can, however, be excluded, since noise in the reconstruction can disturb the particle signature and hence be removed from the estimation. This effect can explain the unpaired particles in Fig. 7. By rejecting false particles, additional calculations will be easier to perform, since false pairing of particles will be avoided.

To compare our results with previously reported results, the numerical aperture and wavelength in the different setups needs to be considered. These parameters govern the axial resolution of the imaging system. The uncertainty will scale with the axial resolution, so a system with a high numerical aperture has a smaller uncertainty. A high numerical aperture accepts a larger set of angles and hence has access to more information from the scattered light. More information creates a more accurate reconstruction. This comparison should not be seen as a benchmark test, but as a rough overview. Compensation for different numerical apertures and wavelengths can be done by scaling. A unitless uncertainty, ${\sigma }_{\text{scaled}}$, is then formulated as

6. CONCLUSIONS

We have showed by simulations and experiments that it is possible to determine the position of a particle by tracking the wavefront curvature from the scattered light. By using off-axis digital holography, the complex amplitude of the light can be recorded and used to reconstruct a volume of the imaged scene with information of both the intensity and phase. From this volume, the wavefront curvature could be estimated by calculating the lateral phase gradients at different axial positions for each particle. The axial position is estimated to be where the phase gradients change sign. From simulations it was found that the estimated position coincides with the true position. Compared to the maximum centroid metric that has a shift, this metric returns the true position without calibration. Experimental measurements were done and it was found that the estimated phase gradients show the same behavior as in the simulations. The phase gradients have a distinct zero-crossing that indicates the true position of the particle. A sample containing the particles was translated 100 μm along the optical axis and holograms were recorded both before and after the translation. Positions from each hologram were estimated and particle pairing between the recordings was performed. The mean displacement distance was found to be 105.0 μm with a standard deviation of 25.3 μm. This was in-line with standard deviations found from simulations. In the simulations, the standard deviation ranged from 2.6 μm to 39.6 μm for noise ratios ranging from 0%to 4% of the maximum bit depth. The phase-gradient method presented in this paper will reject noise from showing up as false particles in the estimation.