1. Introduction

Digital holography is a powerful technique enlarging considerably the depth of investigation of optical microscopy. In one single digital hologram, the different objects can be digitally refocused in the volume without any mechanical scanning as in conventional microscopy. Post-processing of digital holograms is an indispensable step for the analysis of the holographic information. They include among other the extraction of the complex amplitude and the numerical reconstruction [1], processes to correct phase aberrations [2], analysis of living cells [3-5], 3D velocity measurements [6] or furthermore automated 3D object detection and classification [7].

In any automatic object detection process, a major limitation is the number of objects simultaneously under observation. When working with high concentration, the probability that an object goes in front of another increases and it results overlapping of objects. During the recording, those overlapped objects are grabbed and look like an aggregate, considered by processing as one big object instead of the individual ones. An error in the detection/counting process is therefore resulting. In digital holographic microscopy, this problem is recurrent in practical situations when huge quantities of species have to be dynamically analyzed. It is of high importance to overcome this limitation of digital holographic microscopy (DHM) when working with biological samples (like red blood cells or alga cell) for which highly concentrated samples should be considered to get close to realistic situation.

There are two categories of aggregates: those ones constituted by particles in contact (touching particles) and those created by particles lying in different focus planes. In the case of particles in contact, where the constituting objects are approximately in the same focus plane, different separation techniques of images processing are available in the literature. These techniques can be partitioned into model-free and model-based methods.

A well-known model-free approach is the watershed transformation [8], which is a segmentation technique able to detect the ridgelines in an image, by considering its gray-levels as a topographic surface. By combining the watershed segmentation with the distance transform or with the gradient magnitude, clustered nuclei of cells can be divided as demonstrated in several biological studies [9-10]. However, this technique suffers often from over-segmentation due to noise and local irregularities in the image gradient. To overcome this issue, the definition of markers prior to the segmentation is mandatory in many applications. Those markers can be extracted by filtering out all local extrema, irrelevant for the segmentation problem. In [11], the extended h-maxima transform [12] permits to mark cells recorded in in-vitro cultures with inverted light microscopy. More recently, the marker extraction based on the h-minima transform [12] allows segmenting clustered nuclei of cells in fluorescent microscopy images [13] and cervical cells in light microscopy images [14].

Model-based methods can be parametric or non-parametric. An example of parametric method is the Hough transform [15] that can locate well-defined geometry (by analytical curves) in an image through a parameter space. This method has been extended to the generalized Hough transform [16] to detect arbitrary non-analytical shapes by constructing a mapping between the image space and the Hough transform space. In [17], the original circular Hough transform has been modified to detect nuclei of cells in phase contrast microscopy. Non-parametric model-based methods search radial symmetries along the gradient direction. In [18], by iteratively estimating the voting direction (radial or tangential), centroids of nuclear regions can be detected, which permit the partitioning of the regions through Voronoi diagrams. Instead of detecting center of mass of objects, another strategy consists of detecting points of maximum curvature along the contour. In [19], this approach combined with a Delaunay triangulation is used to segment clumps of cell nuclei.

Digital holography offers the possibility to handle aggregates with particles lying in different focus planes. In [20], preliminary results are shown to analyze occluded spherical objects along the optical axis using in-line digital holographic microscopy. By computing the square of the imaginary part of the complex amplitude as a function of the reconstructed distance, local maxima can be slightly highlighted along axial traces revealing the depth positions of the occluded objects. However, the method is only demonstrated for spherical particles of the same sizes, and a robustness estimation of the procedure is not given. Another approach to analyze overlapped objects by digital holography consists to deconvolve the volumetric reconstruction with an optimal kernel that permits to emphasize objects responsible of the scattered patterns [21]. If the sizes of the objects are comparable to the wavelength of light, bright features in the deconvolution reconstruction should lie on the objects themselves. This approach can then only be applied to very small particles, around 1µm. On the other hand, by combining digital holography and Mie scattering theory, recorded holograms of clusters of colloidal spheres can be analyzed to extract their 3D positions. In [22], different clusters of touching colloidal spheres have been fitted to measure translational, rotational and vibrational dynamics. In [23], the same approach has been extended to handle rigid clusters containing up to 6 particles and non-rigid configuration of particles bound to the surface of an emulsion droplet. This method requires initial values for all the model parameters. Despite very good agreement obtained between the recorded hologram and the Mie scattering model, the proposed method is difficult to apply with high concentration of particles [7] due to the high number of parameters and a high computation times. Moreover, as mentioned in [23], the algorithm requires the numbers of spheres in the hologram to be known a priori, which is not the case in practical situations.

In this paper, we investigate the case of aggregates of overlapped particles lying in different focus planes appearing during the recordings of digital holograms. For that purpose, we develop a new technique based on a complete analysis of the evolution of the focus planes that permits to partition the aggregate. Our method allows the handling of aggregates composed of objects that have different shapes, and which lie in different focus planes. The assumption that aggregates are composed of objects having one focus plane is done in this study, which is often the case in practical situations. The algorithm extracts the aggregate border and analyzes the variation of the focus plane along the boundary of the aggregate. Thanks to its simplicity, the developed method can be implemented automatically in any particle detection process and is applicable in many cases. Moreover, it does not require any a priori on the parameter cluster (number, sizes or shapes of the constituting objects) before separation. The robustness of the proposed method is quantified. As shown in the following, aggregates with constituting particles having an overlapping area up to 80% can be processed through the proposed method. This separation procedure can considerably increase the amount of detected particles in high concentrated fluxes.

The paper is organized as follows: in section 2, the developed separation procedure is described through a concrete example. The robustness of the technique is demonstrated by simulations in the first part of section 3 while in the second part, experimental results are exposed to demonstrate the applicability of the developed process in practical situations. Conclusions are given in Section 4.

2. Separation of overlapped particles

The developed method consists to separate aggregates into its constituting particles by computing their positions in depth. The determination of the focus plane of a particle is performed thanks to a refocusing criterion [24-25]. If applied to an aggregate, this criterion will lead to a non-significant value, as there are several focus planes corresponding to the aggregate. To obtain a correct focus plane determination, any application of the refocusing criterion has to be performed on a region containing only one object with a well-defined focus distance. Based on this assumption, the proposed technique consists to generate all around the detected aggregate a set of succeeding regions of interest (ROI) on which the refocusing criterion can be applied. It will result a set of focus planes, one for each ROI, which are merged to provide, by the analysis of the refocus distance distribution, the number of particles constituting the aggregate. The method is described in details in the following section through a concrete example.

2.1 Principle of the developed technique

Consider an aggregate that can be a phase or an amplitude objet. The following analysis can be done in both cases. Figure 1(a) illustrates the amplitude of a simulated aggregate of 4 spherical particles of the same size, with an overlapping area (defined as the area between the overlapped particles) of 64%. With respect to the recorded plane, the particles are defocused respectively by distances of 200, 150 100 and 50µm in a clockwise direction (see Section 3.1 for the description of the simulated aggregate). A 3D view is represented on Fig. 1(b) to highlight the way how the aggregate is constituted.

 

Fig. 1 (a) Intensity image of a simulated aggregate of 4 particles, overlapping each other with an area of 64%. (b) 3D view of the position of the particles inside the simulated volume. The z axis is parallel to the main viewing direction (optical axis).

Download Full Size | PDF

The first post-processing consists of thresholding the object. This is performed by the Otsu method [26] that chooses the threshold level that minimizes the intra-class variance of the black and white pixels. A binary shape covering the object is then obtained with a first computable border. This border extraction can be improved by morphological operations to clean and smooth it. This step is crucial for the outcome of the separation procedure when applied to complex aggregates but can be skipped when working with simple cases. Indeed, the diffraction patterns resulting from the recording of out-of-focus objects create non smoothed complex borders (see Fig. 1(a)). This effect is more sensitive as the overlapping area of particles increases. Therefore, to deal with complex aggregates, smoothing process that is depending on the nature of the extracted border is implemented. This process consists of an opening operation [12]: erosion followed by dilation applied on the binary shape that covers the aggregate. Morphological opening removes completely regions of the object that cannot contain a shape, referred as structuring element (SE), used to probe the object. This element can be of different size and of different type (disk, diamond, square, rectangle, octagon or any arbitrary shape). It has been demonstrated that a generic disk element of size 10 pixels is enough for all the tested aggregates (see next section). Once the object border is smoothed, the border can be extracted. Then each perimeter pixel coordinate is ordered and labeled to form a border segment. The leftmost pixel is chosen as starting point and a clockwise direction is then followed.

As seen on Fig. 1(a), the diffraction rings generated during the propagation of the complex amplitude are radial to the border of the aggregate. In the next section, a theoretical approach is provided to justify the radial propagation of the complex amplitude. Perpendicularly to the border, successive thin rectangular ROIs are generated to cover the complete aggregate, as illustrated for one ROI in blue in Fig. 2(a) . The smoothing procedure by a morphological opening permits the automatic computation of these normal ROIs. Moreover, it avoids strong orientation variations between two adjacent ROIs, which would be the case if the border is noisy and not smooth.

 

Fig. 2 (a) Simulated aggregate of 4 spherical particles. The border in red is extracted from the smoothed border, with an opening operation, of the segmented aggregate. Perpendicularly to this border, successive rectangular ROIs are generated, e.g. in blue, for the first ROI (b) Based on this ROI, a mask is created in which the refocusing criterion is computed.

Download Full Size | PDF

Based on each ROIs, masks are created, as shown in Fig. 2(b), on which the refocusing criterion is computed. It is expected that this criterion will determine correctly the focus planes of the successive ROIs covering the elementary particles. Each ROI is characterized by two parameters: the length and the width. The length is chosen in such a way that the diffraction patterns perpendicular to the border are contained inside the ROI. The width is fixed to recover a small part of the elementary particles. It is also small enough to avoid significant curvature changes when the mask is covering a part of an elementary particle. Tests have been performed to quantify the variation of the computed focus plane as a function of the width of the ROI (for the aggregate of Fig. 2(a)). For very narrow widths (< 5pixels), the diffraction information is not sufficient to reconstruct correctly the focus planes. For large width (> 40pixels), the ROI could recover more than one particle, resulting in a loss of the accuracy of the estimated focus distances. For width between 5 and 30 pixels, the estimated focus distances are acceptable with more accurate results with the smaller width as might be expected. However, the computational time increases with decreasing width. It is then necessary to find some trade-off. In practice, we found that the optimal focus distances are obtained for values of width between 10 and 20 pixels (total image size: 1024x1024 pixels). In Fig. 2(a), the width of the ROI is fixed to 15 and the length to 150 pixels.

Once all the masks are created, the refocusing criterion is computed inside each ROI. As the simulated particles are amplitude objects, the focus plane is reached for the minimum of the criterion [24]. The evolution of the computed focus plane as a function of the pixel perimeter is then plotted, as illustrated on Fig. 3 .

 

Fig. 3 Estimated focus distance as a function of the border segment. The four highlighted horizontal steps correspond to the focus planes of the four simulated spherical particles constituting the aggregate.

Download Full Size | PDF

The four plateaus highlighted on the Fig. 3 correspond to the focus distance of the four simulated spherical particles of the aggregate. Thanks to this representation, the part of the aggregate corresponding to each particle is directly extracted. The computed focus distances between the steps are obtained inside ROIs covering two particles. In those transition regions, the focus distances are thus non-significant and have to be filtered out to determine the number of particles of the under test aggregate. The histogram of the number of border pixels as a function of focus distance can be used to determine the number of particles constituting the aggregate. Their respective focus distances can easily be extracted by a threshold, as shown by Fig. 4 . In our simulated example, we obtained that the four spherical particles are refocused to respectively 200, 149, 99 and 50µm. The separation of the aggregate into its constituting particles is achieved.

 

Fig. 4 Number of perimeter pixels corresponding to each computed focus plane.

Download Full Size | PDF

The mathematical development of the next section shows that the computed ROI are perpendicular to the aggregate border.

2.2 Determination of the orientation of the generated ROIs

In order to determine the suitable shape of the elementary ROIs, we first derive the classical result that the Kirchhoff-Fresnel propagation is equivalent to a Schrödinger-type equation. Consider the complex amplitude distribution $u(x,y,0)$located in a plane$P(z=0)$, propagated up to a distance$d$in a plane$P\text{'}(z=d)$. For that purpose, the Kirchhoff-Fresnel free-space propagation operator (FPO) in the paraxial approximation is used. Thanks to the operator formalism described in [27], it is expressed by:

Considering a small propagation distance$\epsilon$, the quadratic phase operator can be developed until the first order and Eq. (1) becomes:

As a consequence of this Schrödinger equation formulation, the probability conservation equation can be applied to the conserved complex amplitude distribution$u_0(x,y)$:

In the next section, simulations are performed to quantify the robustness of the method. Experiments are exposed to demonstrate its efficiency in practical situations.

3. Results

3.1 Simulations

To quantify the robustness of the method, different types of aggregates were simulated. A cluster generator that can simulate two, three or four particles aggregate was implemented. The aggregate generation procedure can be easily extended to any number of particles. The particles constituting a simulated aggregate can be of various types (spherical, ellipsoidal, cubic or rectangular parallelepiped particles), of different size and orientations. The spherical and ellipsoidal shapes could simulate many real microscopic organisms while the square and rectangular shapes can represent debris that can be present together with organisms in real situations. Each particle is simulated individually in its focus plane. It is then defocused by applying the Fresnel free-space propagation equation [30]. The complex amplitudes of the defocused particles are then recombined to form the aggregate. As explained in the section 2.1, when the overlapping increases, the aggregate becomes more difficult to separate. Therefore, this overlapping area is used as the parameter to quantify the robustness of the method. An aggregate is considered as well separated if, after the application of the separation procedure, each constituting particle can be reconstructed within its depth of field (DOF), defined by$DOF=\lambda /N{A}^2$. The optical system used to perform the experiments, as explained in the next section, uses a laser diode with $\lambda =635nm$and x10 microscope lenses with $NA=0.3$ giving a depth of field of$7\mu m$.

Figure 5 shows a non-exhaustive list of several types of aggregates generated by the cluster generator, with different tested combinations of particle type and orientations. For each aggregate, the overlapping area is increased until the separation process reaches its limit, quantifying thereby the robustness of the developed procedure. In Figs. 5(a)–5(c) are represented aggregates composed of same particle type, respectively spherical, ellipsoidal and cubic ones. In Figs. 5(d)–5(f) aggregates composed with mixed particles are illustrated. For each aggregate, the elementary particles are defocused respectively by 200, 150 100 and 50µm in a clockwise direction. Table 1 gives the reconstructed focus planes as a function of the overlapping area for the aggregate of Figs. 5(a)–5(c). Those aggregates are composed of same types of particles: spherical, ellipsoidal or cubic. The areas between same types of particles composing one aggregate are made identical for each particle. As shown by the Table 1, the method is able to separate the particles of aggregates with overlapping area up to about 80%. Up to this value, the computed focus distances of each particle are within the depth of field, given a maximum error of 4%. This error is defined by $E=|A-B|/A$ where$A$ is the simulated focus distance and$B$is the estimated focus distance obtained with the separation procedure. For aggregates with a larger overlap, the error of the reconstructed focal distance is larger than 7µm and is outside our accepted threshold. By simulations, the robustness of the developed separation process is then established for overlapping area up to about 80% for aggregates composed of identical particles.

 

Fig. 5 Intensity images of different types of aggregates generated by the cluster generator, composed of (a) 4 spherical particles, (b) 4 prolate ellipsoidal particles (2 rotated by a 45° and 2 by a −45° from the x-axis), (c) 4 cubic particles, (d) 1 spherical, 1 ellipsoidal, and 2 cubic particles (whose one is rotated by a 45° from the x-axis), (e) 2 spherical and 2 cubic particles and (f) 2 spherical (of two different sizes) and 2 ellipsoidal particles (whose one is rotated by a 45° from the x-axis).

Download Full Size | PDF

Table 1. Reconstructed focus planes as a function of the area between overlapped particles^a

View Table

For mixed particles aggregates of Figs. 5(d)–5(f), for which the particles have different shapes, the overlapping area is different for each particle. For example, in the simulated aggregate of Fig. 5(d), the overlapping areas are respectively 32.9%, 60.1%, 41.9% and 30.2% for particles in a clockwise direction starting from the upper left particle. Each of those areas changes differently from the others. For mixed particles aggregates of Figs. 5(d)–5(f), each particle is sequentially moved to increase its overlapping area while the others are maintained constant, until the separation method reaches its limit. The robustness is quantified and it has been obtained that a particle of an aggregate is correctly identified if its overlapping area is lower than 75-80%, depending of particles types.

This section demonstrates the capability of the separation procedure to handle different type of simulated aggregates. The robustness of the technique is quantified with respect to the overlapping areas of each particle that can reach up to about 80%.

In the next section, experiments are performed to demonstrate the capacity of the procedure to separate aggregates in practical situations.

3.2 Experiments

The developed technique is applied to real holograms to show the capability of the method in actual situations. Tests were performed with both amplitude and phase objects. A first test is achieved by recording opaque polyethylene spherical particles flowing in a micro-channel. The digital holographic setup used in this experiment is a Mach-Zehnder interferometer in a microscope configuration working with a partial spatial coherent source [7]. The field of view is 720µm x 720µm recorded on a camera with 1024 x 1024 pixels. The x10 microscope lenses have a numerical aperture of NA = 0.3. The diameters of particles are ranging from 75 to 100µm. The micro-fluidic experimental device consists of a straight micro-channel µ-Slide$I^{0.8}$Ibidi® (height: 800µm). To perform the flow measurements, we use a programmable KDS Legato® 270P push/pull pump with a user defined program configuration, allowing a wide flow range.

For opaque particles, the detection of the aggregate is performed on the intensity image extracted from the recorded hologram. It has been shown that the focus plane of an amplitude object is reached for the minimum of a refocusing criterion defined in [24-25]. Figure 6 shows three different aggregates and the corresponding separation into their constituting particles with the separation method. By analyzing the evolution of the focus planes around the aggregate border, different levels are highlighted as shown in Fig. 3, where the focus planes of the different particles can be directly extracted. Figures 6(a), 6(d) and 6(g) show the intensity image of aggregates of respectively two, three and four opaque spherical particles of different sizes. Thanks to the separation method, the focus planes of each particle of the aggregates are computed as shown by Figs. 6(b)–6(c) for the first aggregate, by Figs. 6(e)–6(f) for the second one and by Figs. 6(h)–6(i) for the third one. For real data, as the focus planes are not known a priori, the refocus distance errors on each particle are estimated by computing the standard deviation of the mean recomputed focus planes of all the ROIs recovering the particle. The obtained errors are below our acceptable threshold.

 

Fig. 6 Separation method applied on aggregates constituted by opaque polyethylene spherical particles of different sizes. (a) Intensity image of an aggregate of two particles in the recorded plane (b) Automated separation and refocusing of the bottom particle at −370µm (error: 5µm). (c) Automated separation and refocusing of the top particle at 320µm (error: 6µm). (d) Intensity image of an aggregate of three particles in the recorded plane (e) Automated separation and refocusing of the top particle at −135µm (error: 3µm). (f) Automated separation and refocusing of the bottom particle at 95µm (error: 2µm). (g) Intensity image of an aggregate of four particles in the recorded plane (h) Automated separation and refocusing of the top particle at −125µm (error: 2µm). (i) Automated separation and refocusing of the bottom particle at −265µm (error: 4µm).

Download Full Size | PDF

Other tests were performed with red blood cells (RBCs). This problem is of high interest as blood is made of about 50% of RBCs, which means that the volume fraction of RBC in real blood is very high. For this reason, blood must be strongly diluted to access information of individual RBC (typically 0.1% in volume fraction). Even with high dilution, the probability to have an aggregate of RBC is high. Those clusters corrupt the determination of the 3D locations of RBCs. The application of our separation method on RBCs cluster is a first step towards highest concentration monitoring. RBCs are oval and biconcave disks with diameters ranging from 6 to 8µm. The thickness at the thickest location is approximately 2-2.5µm while the thickness in the center is about 1µm.

The optical setup has been set to yield a field of view of 215 µm x 215 µm obtained with a 40x microscope lenses, with NA = 0.6. The resulting DOF is then equal to 1.8µm for this setup. A suspension of RBC is injected in a shear-flow chamber made of two glass plates spaced by 170µm. The bottom disk is fixed while the top one is rotating, creating a linear velocity profile of the fluid inside the channel [31]. Holograms are acquired at a frame rate of 24 holograms per second. The main objective of those experiments is to get the depth position of each RBC and the evolution in time of the RBC suspension.

Figures 7(a) –7(b) represent respectively the intensity and the phase image of an aggregate of two RBCs overlapping each other. The detection and segmentation process, as explained in Section 2.1, is performed on the phase image and is illustrated in Fig. 7(c). We can see here that an automatic detection process without separation will consider the aggregate as one big object.

 

Fig. 7 Cluster of two RBCs in the recorded plane (a) Intensity image (b) Phase image (c) Segmentation of the phase image of the aggregate showing that one big object is detected.

Download Full Size | PDF

Thanks to the separation method, the RBC on the left is automatically separated from the aggregate and refocused at −24µm of the recorded plane, as illustrated by Figs. 8(a) –8(b), while the RBC on the right is refocused at −76µm from the recorded plane, as seen on Figs. 8(c)–8(d). The computed errors are in our accepted margin of 1.8 µm. Once the cluster is separated, morphological features of RBCs can be extracted in their focus plane to study their 3D orientation in the channel [32-33].

 

Fig. 8 Separation procedure applied to RBCs cluster (a) Automated separation and refocusing of the left particle at −24µm (error: 0.8µm) - Intensity image (b) Automated separation and refocusing of the left particle at −24µm (error: 0.8µm) - Phase image (c) Automated separation and refocusing of the left particle at −76µm (error: 1.5µm) - Intensity image (d) Automated separation and refocusing of the left particle at −76µm (error: 1.5µm) - Phase image.

Download Full Size | PDF

The difficulty of working with RBCs lies in the fact that depending on their orientation, their projected shape might be non-symmetrical. Different clusters shapes are then possible during the hologram recording as illustrated in Fig. 9 . On this cluster, the two RBCs have two different sizes and share a higher overlapping area than in the previous example as observed on the intensity image and phase image of Figs. 9(a)–9(b). After applying our separation procedure, the RBC on the right is automatically separated from the cluster and refocused at −90µm from the recorded plane, as illustrated by Figs. 9(c)–9(d), while the RBC on the left is refocused at −42µm from the recorded plane, as seen on Figs. 9(e)–9(f).

 

Fig. 9 Cluster of two RBCs in two different orientations (a) Intensity image of the cluster in the recorded plane (b) Phase image of the cluster in the recorded plane (c) RBC on the right separated and refocused at −90µm (error: 2.4µm) – Intensity image (d) RBC on the right separated and refocused at −90µm (error: 2.4µm) – Phase image (e) RBC on the left separated and refocused at −42µm (error: 1.2µm) – Intensity image (f) RBC on the left separated and refocused at −42µm (error: 1.2µm) – Phase image.

Download Full Size | PDF

For this cluster, the right cell is deformed and elongated along the flow direction (from bottom to top on the image). Furthermore, this cell lies in a plane parallel to the optical axis and exhibits a prolate ellipsoidal shape, while the other RBC lies in a plane perpendicular to the optical axis. The RBC on the right has then a higher optical thickness than the left one. It is the reason why its error in the refocused z position is higher than the left one. Even in this complex case where the obtained error is slightly above our limit, the estimated z position is acceptable for morphological feature extraction as we can see on Fig. 9(d).

These experiments confirm the capability of the developed separation method in practical situations.

3. Conclusions

In this article, a new method to separate aggregate of particles is proposed in digital holographic microscopy. The developed procedure is based on an accurate analysis of the evolution of the focus plane around the aggregate’s border. Aggregates of amplitude objects as well as phase object can be separated with the proposed technique, which can moreover deal with many types of aggregate shape. Numerical simulations permit to quantify the robustness of the method. It has been demonstrated that the overlapping area can be increased up to a level of around 80% recovering the correct focus planes of the objects inside their depth of field. Experiments were performed and demonstrated that the method can be applied in practical situations.