Method for counting particles in a sample by means of lensless imaging
11199488 · 2021-12-14
Assignee
- Commissariat A L'energie Atomique Et Aux Energies Alternatives (Paris, FR)
- Horiba Abx Sas (Montpellier, FR)
Inventors
- Pierre Blandin (Coublevie, FR)
- Anais Ali-Cherif (Clermont Ferrand, FR)
- Estelle Gremion (La Tour du Pin, FR)
- Sebastien Raimbault (Argelliers, FR)
- Olivier Cioni (Grenoble, FR)
- Aurelien Daynes (Montpellier, FR)
Cpc classification
G03H1/0866
PHYSICS
G01N2015/1454
PHYSICS
G03H1/0443
PHYSICS
International classification
Abstract
The invention relates to a method for counting particles, particularly blood cells, in a sample, using a lensless optical imaging device. The sample is arranged between a light source and an image sensor. The sample is illuminated by a light source and an image is acquired by the image sensor, said image sensor being exposed to a light wave called an exposition wave. A digital propagation operator is applied to the acquired image so as to obtain a complex amplitude of the exposition wave according to a surface facing the image sensor. An image, called a reconstructed image, is formed from the modulus and/or the phase of said complex amplitude, on which image the particles to be counted appear in the form of regions of interest. The method then comprises a step of selecting the regions of interest corresponding to the particles to be counted.
Claims
1. A method for counting particles, known as particles of interest, arranged in a sample, the method comprising the following steps: a) illumination of the sample with the aid of a light source, the light source emitting an incident light wave propagating toward the sample; b) acquisition, with the aid of an image sensor, of an image (I.sub.0) of the sample, formed in a detection plane (P.sub.0), the sample being arranged between the light source and the image sensor, each image being representative of a light wave, known as the exposition wave, to which the image sensor is exposed under the effect of said illumination; c) based on the image acquired (I.sub.0) during step b), application of a propagation operator (h) so as to calculate a complex expression (A.sub.z) of the exposition wave, along a reconstruction surface (P.sub.z) extending opposite the detection plane (P.sub.0); d) formation of an image (I.sub.z, M.sub.z, φ.sub.z), known as the reconstructed image, representative of a distribution of the modulus or of the phase of the complex expression along the reconstruction surface; e) segmentation of the image formed during step d) to obtain a segmentation image (I.sub.z*), containing regions of interest (ROI) spaced apart from each other, all or some of the regions of interest being associated with at least one particle of interest; f) determination of a size of regions of interest, and classification of said regions of interest as a function of their size according to at least one of the following classes: a class such that the region of interest contains a single particle of interest; a class such that the region of interest contains a whole number, greater than 1, of particles of interest; and g) counting of the particles of interest based on the regions of interest which underwent a classification during step f).
2. The method as claimed in claim 1, in which step f) comprises, prior to the classification: a determination of at least one morphological criterion for each region of interest; a selection of the regions of interest (ROI) as a function of the morphological criterion, the classification being done for the selected regions of interest.
3. The method as claimed in claim 2, wherein at least one morphological criterion is chosen from among a diameter, a size, or a form factor of each region of interest.
4. The method as claimed in claim 1, wherein step f) comprises taking into account a reference size (S.sub.ref), representative of a predetermined number of particles of interest, the classification being done as a function of said reference size.
5. The method as claimed in claim 1, wherein step f) comprises the following sub-steps: fi) a first classification of the regions of interest, performed as a function of a first reference size (S.sub.ref−1), the first reference size corresponding to a predetermined number of particles of interest; fii) after the first classification, a determination of a second reference size (S.sub.ref−2), based on the regions of interest classified, during sub-step fi), as containing a predefined number of particles of interest; and fiii) a second classification of the regions of interest, performed as a function of the second reference size (S.sub.ref−2) determined during sub-step fii).
6. The method as claimed in claim 1, wherein step f) comprises, prior to the classification: a calculation of a signal to noise ratio (S/B) for several regions of interest (ROI); and a selection of the regions of interest (ROI) as a function of the signal to noise ratio (S/B) calculated for each of them, the classification being done for the regions of interest so selected.
7. The method as claimed in claim 1, wherein the sample is mixed with a sphering reagent prior to step b), in order to modify the shape of the particles of interest so as to make it spherical.
8. The method as claimed in claim 1, wherein step c) comprises the following sub-steps: ci) definition of an initial image of the sample (A.sub.0.sup.k=0) in the detection plane, based on the image (I.sub.0) acquired by the image sensor; cii) determination of a complex image of the sample (A.sub.z.sup.k) in a reconstruction surface (P.sub.z) by applying a propagation operator to the initial image of the sample (A.sub.0.sup.k=1) defined during sub-step ci) or to the image of the sample (A.sub.0.sup.k−1), in the detection plane, resulting from the preceding iteration (k−1); ciii) calculation of a noise indicator (∈.sup.k) based on the complex image (A.sub.z.sup.k) determined during sub-step cii), this noise indicator depending on a reconstruction noise affecting said complex image (A.sub.z.sup.k); civ) updating of the image of the sample (A.sub.0.sup.k) in the detection plane (P.sub.0) by an adjustment of phase values (φ.sub.0.sup.k(x,y)) of the pixels of said image, the adjustment being done as a function of a variation of the indicator calculated during sub-step ciii) according to said phase values; and cv) reiteration of sub-steps cii) to civ) until achieving a criterion of convergence, so as to obtain a complex image of the sample in the detection plane (P.sub.0) and in the reconstruction surface (P.sub.z).
9. The method as claimed in claim 1, wherein no image-forming optic is placed between the sample and the image sensor.
10. The method as claimed in claim 1, wherein the particles of interest are blood cells.
11. The method as claimed in claim 10, wherein during step d): the particles of interest are red blood cells or white blood cells, the reconstructed image being representative of a distribution of the modulus of the complex expression (A.sub.z) along the reconstruction surface (P.sub.z); and the particles of interest are platelets, the reconstructed image being representative of a distribution of a phase of the complex expression (A.sub.z) along the reconstruction surface (P.sub.z).
12. The method as claimed in claim 1, wherein step c) comprises an application of a digital focusing algorithm so as to determine a focusing distance corresponding to the distance of reconstruction to which the reconstruction surface (P.sub.z) extends.
13. The method as claimed in claim 12, wherein during step c) the image acquired during step b) is partitioned into different zones, distinct from each other, and in which the digital focusing algorithm is applied in each zone so as to obtain, for each of them, a focusing distance (d.sub.p), step c) also comprising: taking into account a maximum dispersion indicator (δ.sub.max); calculation of a dispersion indicator (δ(d.sub.p)) representing a dispersion of the obtained focusing distances (d.sub.p); a comparison between the calculated dispersion indicator (δ(d.sub.p)) and the maximum dispersion indicator (δ.sub.max); as a function of the comparison: calculation of a mean or a median (
14. The method as claimed in claim 13, wherein the dispersion indicator (δ(d.sub.p)) is a difference between a maximum focusing distance (d.sub.p−max) and a minimum focusing distance (d.sub.p−min), or a variance or a standard deviation of the obtained focusing distances (d.sub.p).
15. The method as claimed in claim 12, wherein during step c) the image acquired during step b) is partitioned into different zones, distinct from each other, and in which the digital focusing algorithm is applied in each zone so as to obtain, for each of them, a focusing distance (d.sub.p), step c) also comprising: taking into account a lower focusing limit (d.sub.inf); taking into account an upper focusing limit (d.sub.sup); a comparison of each focusing distance (d.sub.p) with the lower focusing limit (d.sub.inf) and with the upper focusing limit (d.sub.sup); as a function of the comparison: if the focusing distance is between the lower focusing limit (d.sub.inf) and the upper focusing limit (d.sub.sup), calculation of a mean or a median (
16. The method as claimed in claim 1, wherein step e) comprises: a determination of a mode of the histogram of the reconstructed image (I.sub.z, M.sub.z, φ.sub.z), corresponding to an intensity of a maximum number of pixels of the reconstructed image; implementing an Otsu thresholding in order to obtain a segmentation threshold based on which the segmentation of the reconstructed image is done; and a comparison between the mode of the histogram and the segmentation threshold, such that when the segmentation threshold is greater than 90% of the mode of the histogram the segmentation threshold is set at a percentage between 80% and 85% of the mode of the histogram.
17. The method as claimed in claim 1, wherein step e) comprises a determination of a criterion of quality (Q) of the segmented image (I.sub.z*), comprising: taking into account a predetermined maximum area (A.sub.max); a determination of the area (A) of all of the regions of interest (ROI) of the segmented image (I.sub.z*), a comparison of the area (A) of all of the regions of interest with the maximum area (A.sub.max); and as a function of the comparison, rejection or validation of the segmented image (I.sub.z*).
18. The method as claimed in claim 17, wherein, if the segmented image (I.sub.z*) is validated, and the area (A) of the regions of interest is situated close to the maximum area (A.sub.max), the method comprises a determination of an auxiliary criterion of quality (Q′), comprising: taking into account a predetermined minimum signal to noise ratio (SNR.sub.min); calculation of a signal to noise ratio (SNR(I.sub.z*) of the segmented image (I.sub.z*), in the form of a ratio between a mean or median intensity (I.sub.ROI) in the regions of interest (ROI) and a noise level (N.sub.out) of the segmented image outside the regions of interest; a comparison between the signal to noise ratio of the segmented image (SNR(I.sub.z*)) and the minimum signal to noise ratio (SNR.sub.min); and as a function of the comparison, rejection or validation of the segmented image (I.sub.z*).
19. The method as claimed in claim 1, wherein during step f) the regions of interest whose size exceeds a maximum predetermined size (S.sub.max) are rejected, the method comprising: a determination of the number of regions of interest rejected during step f); taking into account a predetermined level of acceptance (L.sub.max); and a comparison between the number of regions of interest rejected during step f) and the level of acceptance (L.sub.max) as a function of which the counting of the particles of interest, performed during step g), is validated or invalidated.
20. A device for the counting of particles of interest arranged in a sample, the device comprising: a light source to emit an incident light wave propagating toward the sample; a support designed to hold the sample between the light source and an image sensor; and a processor configured to receive an image of the sample acquired by the image sensor and to carry out steps c) to g) of the method as claimed in claim 1.
Description
FIGURES
(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
(9)
(10)
PRESENTATION OF PARTICULAR EMBODIMENTS
(11)
(12) The sample 10 is a sample containing particles, among which are particles of interest 10a which one wishes to count. By particles is meant cells, for example, and in particular blood cells, but they may also be microorganisms, viruses, spores, or microbeads, usually employed in biological applications, or even microalgae. They may also be droplets which are insoluble in a liquid medium 10m, such as oil droplets dispersed in an aqueous phase. Preferably, the particles 10a have a diameter, or are inscribed in a diameter, less than 100 μm, and preferably less than 50 μm or 20 μm. Preferably, the particles have a diameter, or are inscribed in a diameter, greater than 500 nm or 1 μm.
(13) In the example represented in
(14) The sample 10 in this example is contained in a fluidic chamber 15. The fluidic chamber 15 is for example a fluidic chamber of Countess® type with a thickness e=100 μm. The thickness e of the sample 10, along the axis of propagation Z, typically varies between 10 μm and 1 cm, and is preferably between 20 μm and 500 μm. The sample 10 extends along a plane P.sub.10, known as the plane of the sample, preferably perpendicular to the axis of propagation Z. It is maintained on a support 10s at a distance from an image sensor 16.
(15) The distance D between the light source 11 and the sample 10 is preferably greater than 1 cm. This is preferably between 2 and 30 cm. Advantageously, the light source, as seen by the sample, is considered to be pointlike. This means that its diameter (or its diagonal) is preferably less than one tenth, or better one hundredth, of the distance between the sample and the light source. The light source 11 may be a light-emitting diode as shown in
(16) The device may comprise a diffuser 17, arranged between the light source 11 and the diaphragm 18. The use of such a diffuser makes it possible to overcome constraints of centering of the light source 11 with respect to the aperture of the diaphragm 18. The function of such a diffuser is to distribute the light beam produced by the light source along a cone with angle α. Preferably, the angle of diffusion a varies between 10° and 80°. The presence of such a diffuser helps make the device more tolerant of off-centering of the light source with respect to the diaphragm, and homogenizes the illumination of the sample. The diaphragm is not necessary, especially when the light source is sufficiently pointlike, particularly in the case of a laser source.
(17) The light source 11 may be a laser source, such as a laser diode. In this case, it is not useful to associate it with a spatial filter 18 or a diffuser 17. Such a configuration is shown in
(18) Preferably, the spectral emission band Δλ of the incident light wave 12 has a width less than 100 nm. By spectral band width is meant a width at mid-height of said spectral band.
(19) The sample 10 is arranged between the light source 11 and the previously mentioned image sensor 16. The latter extends preferably in parallel or substantially in parallel with the plane P.sub.10 of extension of the sample. The term substantially parallel means that the two elements might not be strictly parallel, an angular tolerance of a few degrees, less than 20° or 10°, being allowed. The image sensor 16 is able to form an image I.sub.0 along a detection plane P.sub.0. In the example shown, it is an image sensor comprising an array of pixels, of CCD type, or a CMOS. The detection plane P.sub.0 preferably extends perpendicular to the axis of propagation Z of the incident light wave 12. The distance d between the sample 10 and the array of pixels of the image sensor 16 is advantageously between 50 μm and 2 cm, preferably between 100 μm and 2 mm.
(20) One will note the absence of any image-forming optic, such as a magnification optic, between the image sensor 16 and the sample 10. This does not preclude the possible presence of focusing micro-lenses in the area of each pixel of the image sensor 16, the latter not having the function of magnification of the image acquired by the image sensor.
(21) Under the effect of the incident light wave 12, the particles present in the sample may create a diffracted wave 13, able to produce interference in the area of the detection plane P.sub.0 with a portion 12′ of the incident light wave 12 transmitted by the sample. Moreover, the sample may absorb a portion of the incident light wave 12. Thus, the light wave 14, known as the exposition wave, which is transmitted by the sample 10 and to which the image sensor 16 is exposed, may comprise: a component 13 resulting from the diffraction of the incident light wave 12 by each particle of the sample; a component 12′ resulting from the transmission of the incident light wave 12 by the sample.
(22) These components produce interference in the detection plane. Thus, the image I.sub.0 acquired by the image sensor 16 comprises interference patterns (or diffraction patterns), each interference pattern being generated by a particle 10a of the sample 10.
(23) A processor 20, such as a microprocessor, is designed to process each image I.sub.0 acquired by the image sensor 16. In particular, the processor is a microprocessor connected to a programmable memory 22 in which a sequence of instructions is memorized for carrying out the image processing and calculation operations described in this description. The processor may be coupled to a screen 24 enabling the display of images acquired by the image sensor 16 or calculated by the processor 20.
(24) As mentioned in connection with the prior art, an image I.sub.0 acquired at the image sensor 16, also known as a hologram, does not allow one to obtain a sufficiently precise representation of the sample observed. One may apply, to each image acquired by the image sensor, a propagation operator h, in order to calculate a quantity representative of the exposition wave 14. Such a procedure, known as holographic reconstruction, in particular allows the reconstructing of an image of the modulus or of the phase of this light wave 14 in a reconstruction plane parallel to the detection plane P.sub.0, and in particular in the plane P.sub.10 of extension of the sample. For this, one makes a convolution product of the image I.sub.0 acquired by the image sensor 16 with a propagation operator h. It is then possible to calculate a complex expression A of the light wave 14 at every point of spatial coordinates (x,y,z), and in particular along a reconstruction surface extending opposite the image sensor. The reconstruction surface may be in particular a reconstruction plane P.sub.z situated at a distance |z| from the image sensor 16, and in particular the plane P.sub.10 of the sample, with: A(x, y, z)=I.sub.0(x, y, z)*h, * denoting the operator produced from the convolution.
(25) In the rest of this description, the coordinates (x, y) denote a radial position in a plane perpendicular to the axis of propagation Z. The coordinate z denotes a coordinate along the axis of propagation Z.
(26) The complex expression A is a complex quantity whose argument and modulus are respectively representative of the phase and the intensity of the light wave 14 to which the image sensor 16 is exposed. The convolution product of the image I.sub.0 with the propagation operator h makes it possible to obtain a complex image A.sub.z representing a spatial distribution of the complex expression A in the reconstruction plane P.sub.z, extending to a coordinate z of the detection plane P.sub.0. In this example, the detection plane P.sub.0 has for its equation z=0. The complex image A.sub.z corresponds to a complex image of the sample in the reconstruction plane P.sub.z. It likewise represents a two-dimensional spatial distribution of the optical properties of the exposition wave 14. The propagation operator h has the function of describing the propagation of the light between the image sensor 16 and a point of coordinates (x, y, z), situated at a distance |z| from the image sensor. It is thus possible to determine the modulus M(x,y,z) and/or the phase φ(x, y, z) of the exposition wave 14, at this distance |z|, known as the distance of reconstruction, with:
M(x,y,z)=abs[A(x,y,z)] (1);
φ(x,y,z)=arg [A(x,y,z)] (2).
(27) The operators abs and arg denote respectively the modulus and the argument.
(28) In other words, the complex expression A of the exposition wave 14 at every point of the spatial coordinates (x, y, z) is such that: A(x, y, z)=M(x, y, z)e.sup.jφ(x, y, z) (3). It is possible to form images M.sub.z and φ.sub.z representing respectively a distribution of the modulus or of the phase of the complex expression A in a surface extending opposite the detection plane P.sub.0. Such a surface may be in particular a plane P.sup.z situated at a distance |z| from the detection plane P.sub.0, with M.sub.z=abs(A.sub.z) and φ.sub.z=arg(A.sub.z). The previously mentioned surface is not necessarily planar, but it may extend substantially parallel to the detection plane and is preferably a plane P.sup.z parallel to the detection plane. In the rest of the description, the image obtained from the modulus and/or the phase of the complex image A.sub.z is called the “reconstructed image” and denoted by I.sub.z.
(29) The inventors have developed a method for counting particles 10a present in the sample by a procedure described in connection with
(30) Step 100: illumination of the sample. During this step, the sample is illuminated by the light source 11.
(31) Step 110: acquisition of an image I.sub.0 of the sample 10 by the image sensor 16, this image forming a hologram. One interesting feature of the lensless configuration, as represented in
(32) Step 120: calculation of a complex image A.sub.z in a reconstruction plane P.sub.z. The complex image carries information about the phase and amplitude of the exposition wave 14 to which the image sensor 16 is exposed. The reconstruction plane is the plane P.sub.10 of extension of the sample. Step 120 may be carried out by applying the propagation operator h, previously described, to an image resulting from the acquired image I.sub.0. However, the application of the propagation operator to the acquired image may result in a complex image A.sub.z containing substantial reconstruction noise, often known as a twin image. In order to obtain a usable complex image, limiting the reconstruction noise, iterative algorithms may be implemented. One of these algorithms is described below, in connection with
(33) From the complex image A.sub.z, one may obtain a reconstructed image I.sub.z based on the modulus M.sub.z and/or the phase φ.sub.z of the exposition wave 14, in the reconstruction plane P.sub.z. The reconstruction distance with respect to the image sensor 16 is determined either a priori, knowing the position of the sample 10 with respect to the image sensor 16, or on the basis of a digital focusing algorithm, whereby the optimal reconstruction distance is the one for which the reconstructed image I.sub.z is the most distinct. Digital focusing algorithms are known to the person skilled in the art.
(34)
(35) Step 130: segmentation. The image formed during step 120 is subjected to a segmentation, in order to isolate the regions of interest ROI corresponding to particles. By image segmentation is meant a partitioning of the image so as to regroup the pixels, for example as a function of their intensity. The image segmentation results in a segmentation image I.sub.z* in which the regions of interest ROI, spaced apart from each other, are delimited, each region of interest possibly corresponding to a particle, or to a cluster of particles, as described below. Different methods of segmentation are known to the person skilled in the art. For example, one may apply an Otsu thresholding, consisting in determining a value of an intensity threshold from the histogram of the image, the threshold allowing an optimal separation of the pixels according to two classes: one class of pixels representing the regions of interest and one class of pixels representing the background of the image.
(36) In this example, the segmentation image I.sub.z* obtained in this step results from a segmentation of the image M.sub.z of the modulus of the complex amplitude A.sub.z of the exposition wave 14. As a variant, the image obtained in this step results from a segmentation of the image φ.sub.z of the phase of the exposition wave. The segmentation image I.sub.z* may likewise combine the regions of interest appearing after the respective segmentations of the image M.sub.z of the modulus as well as the image φ.sub.z of the phase of the exposition wave.
(37) A counting of the particles of interest by a simple counting of the regions of interest appearing in the segmentation image I.sub.z* may provide a first order of magnitude as to the quantity of particles of interest. But the inventors have discovered that such a counting would yield rather imprecise results. In fact, certain regions of interest correspond to particles different from the particles that one wishes to count. Moreover, certain regions of interest bring together several particles of interest, which form clusters. This is due to the fact that certain particles of interest are close to each other, and are merged into the same region of interest. In
(38) Step 140: selection. During this step, a selection is done for the regions of interest that were determined during the previous step as a function of their morphology, that is, as a function of a morphological criterion, which might be the area, the shape, or the size. This makes it possible to select the regions of interest corresponding to the particles of interest that one wishes to count.
(39) This step may involve a filtering of each region of interest as a function of the previously mentioned morphological criterion, in order to select regions of interest representative of the particles of interest. For example, when the particles of interest are red blood cells, the filtering makes it possible to select the regions of interest having a diameter less than 100 μm, or inscribed in such a diameter. This allows not considering large dust motes or traces, while preserving clusters of red blood cells.
(40) In order to take account of effects of agglutination between particles of interest, a classification of the selected regions of interest may be done as a function of their area and/or their shape, especially based on a reference size S.sub.ref, such as a reference area, the latter being for example representative of a region of interest corresponding to a single particle of interest (singlet). The reference area may be predetermined or obtained from the segmentation image I.sub.z*, for example by calculating a mean, or a median, of the area of each region of interest ROI, and assuming that the majority of the regions of interest correspond to a single particle. The size of each region of interest is then compared to the reference size S.sub.ref, as a function of which each region of interest is assigned a number of particles. By size is meant a diameter or an area or another characteristic dimension. For example, in connection with
(41) The notation ROI.sub.n corresponds to a region of interest which is considered, after the classification step, to comprise n particle(s), n being an integer which is strictly positive.
(42) In an alternative or supplementary manner, prior to the classification, a selection may be done as a function of the form of each region of interest. For example, one determines for each region of interest a largest diameter and a smallest diameter. The ratio between the largest diameter and the smallest diameter makes it possible to quantify a form factor of the considered region of interest. By this selection, one may identify the regions of interest whose shape cannot be considered to be spherical, these regions of interest not being considered as being representative of particles of interest and not being taken into account in the counting. Thus, the selection may involve a comparison between the shape of each region of interest and a shape representative of the particle of interest, or cluster of particles of interest, being detected.
(43) For each region of interest present in the segmentation image I.sub.z*, a signal to noise ratio S/B can be established. Such a signal to noise ratio may be calculated by producing a ratio between a mean value of pixels, known as central pixels, situated at the center of the region of interest, for example 9 central pixels, and the standard deviation calculated for a background of the image, the background of the image corresponding to the reconstructed image without the regions of interest.
(44) According to one embodiment, step 140 may involve the succession of the following operations, described in connection with
(45) This algorithm allows a progressive adjustment of the reference size used as the basis for the classification of the regions of interest. This allows one to obtain a more precise classification.
(46) Thus, step 140 makes it possible to select the regions of interest which are representative of the particles of interest, and/or to assign a number of particle(s) to each region of interest considered as being representative of particles of interest.
(47) Step 150: counting. During this step, one counts the regions of interest selected during step 140, that is, those considered as being representative of particles, taking into account a number of particles possibly assigned to each of them.
(48) Step 120 described above may be realized by producing a convolution between a propagation operator and the image acquired by the image sensor. However, the application of such an operator may give rise to significant reconstruction noise. In order to optimize the reconstruction by limiting the reconstruction noise, step 120 may be carried out by an algorithm such as is described in the patent application FR1652500 filed 23 Mar. 2016. We shall now describe, in connection with
(49) Step 121: propagation of the detection plane toward the reconstruction plane
(50) During this step, the image formed in the detection plane P.sub.0 is available. During the first iteration, an initial image A.sub.0.sup.k=0 is determined from the image I.sub.0 acquired by the image sensor. The modulus M.sub.0.sup.k=0 of the initial image A.sub.0.sup.k=0 may be obtained by applying the square root operator to the image I.sub.0 acquired by the image sensor, in which case M.sub.0.sup.k=0=√{square root over (I.sub.0)}. An arbitrary value, such as 0, is assigned to the phase φ.sub.0.sup.k=0 of the initial image. During the following iterations, the image in the detection plane is the complex image A.sub.0.sup.k−0 resulting from the previous iteration. The image formed in the detection plane P.sub.0 is propagated in the reconstruction plane P.sub.z, by the application of a propagation operator h as previously described, so as to obtain a complex image A.sub.z.sup.k, representative of the sample 10, in the reconstruction plane P.sub.z. The propagation is realized by convolution of the image A.sub.0.sup.k−0 with the propagation operator h.sub.−z, such that:
A.sub.z.sup.k=A.sub.0.sup.k−0*h.sub.−z,
(51) The index −z represents the fact that the propagation is realized in a direction opposite to the axis of propagation Z. One speaks of retropropagation.
(52) Step 122: Calculation of an indicator in several pixels
(53) During this step, one calculates a quantity ∈.sup.k(x, y) associated with each pixel of a plurality of pixels (x, y) of the complex image A.sub.z.sup.k, and preferably in each of its pixels. This quantity depends on the value A.sub.z.sup.k(x, y) of the image A.sub.z.sup.k, or its modulus, at the pixel (x, y) for which it is calculated. It may likewise depend on a dimensional derivative of the image in this pixel, for example the modulus of a dimensional derivative of this image. In this example, the quantity ε.sup.k(x, y) associated with each pixel is a modulus of a difference of the image A.sub.z.sup.k, in each pixel, and the value 1. Such a quantity may be obtained by the expression:
ε.sup.k(x,y)=√{square root over ((A.sub.z.sup.k(x,y)−1)(A.sub.z.sup.k(x,y)−1)*)}=|A.sub.z.sup.k(x,y)−1|
(54) Step 123: establishing of a noise indicator associated with the image A.sub.z.sup.k.
(55) During step 122, one calculates quantities ε.sup.k(x, y) in several pixels of the complex image A.sub.z.sup.k. These quantities may form a vector E.sup.k, whose terms are the quantities ε.sup.k(x, y) associated with each pixel (x, y). In step 123, one calculates an indicator, known as a noise indicator, from a norm of the vector E.sup.k. The quantity ε.sup.k(x, y) calculated from the complex image A.sub.z.sup.k, at each pixel (x, y) of the latter, is summed in order to constitute a noise indicator ε.sup.k associated with the complex image A.sub.ref.sup.k.
(56) Thus, ε.sup.k=Σ.sub.(x, y)ε.sup.k(x, y)
(57) An important aspect of this step is to determine, in the detection plane P.sub.0, phase values φ.sub.0.sup.k(x, y) of each pixel of the image A.sub.0.sup.k in the plane of the sample, making it possible to obtain, during a following iteration, a reconstructed image A.sub.z.sup.k+1 whose noise indicator ε.sup.k+1 is less than the noise indicator ε.sup.k.
(58) During the first iteration, the only information available pertains to the intensity of the exposition wave 14, but not to its phase. The first reconstructed image A.sub.z.sup.k=1 in the reconstruction plane P.sub.z is thus assigned a substantial reconstruction noise, due to the absence of information as to the phase of the light wave 14 in the detection plane P.sub.0. Consequently, the indicator ε.sup.k=1 is elevated. During subsequent iterations, the algorithm carries out a progressive adjustment of the phase φ.sub.0.sup.k(x, y) in the detection plane P.sub.0, so as to progressively minimize the indicator ε.sup.k.
(59) The image A.sup.0.sub.k in the detection plane is representative of the light wave 14 in the detection plane P.sub.0, both from the standpoint of its intensity and its phase. Steps 120 to 125 are designed to establish, in iterative fashion, the value of the phase φ.sub.0.sup.k(x, y) of each pixel of the image A.sub.z.sup.k, minimizing the indicator ε.sup.k, the latter having been obtained for the image A.sub.z.sup.k obtained by propagation of the image A.sub.0.sup.k−1 in the reconstruction plane P.sub.z.
(60) The minimization algorithm may be a gradient descent algorithm or a conjugated gradient descent algorithm, the latter being described below.
(61) Step 124: Adjustment of the value of the phase in the detection plane.
(62) Step 124 is designed to determine a value of the phase φ.sub.0.sup.k(x, y) of each pixel of the complex image A.sub.0.sup.k so as to minimize the indicator ε.sup.k+1 resulting from a propagation of the complex image A.sub.0.sup.k in the reconstruction plane P.sup.z, during the following iteration k+1. For this, a phase vector φ.sub.0.sup.k is established, each term of which is the phase φ.sub.0.sup.k(x, y) of a pixel (x, y) of the complex image A.sub.0.sup.k. The dimension of this vector is (N.sub.pix, 1), where N.sub.pix denotes the number of pixels considered. This vector is updated during each iteration, by the following updating expression:
φ.sub.0.sup.k(x,y)=φ.sub.0.sup.k−1(x,y)+α.sup.kp.sup.k(x,y) where: α.sup.k is an integer, denoted as the “step”, and representing a distance; p.sup.k is a direction vector, of dimension (N.sub.pix, 1), each term p(x, y) of which forms a direction of the gradient ∇ε.sup.k of the indicator ε.sup.k.
(63) This equation may be expressed in vectorial form as follows:
φ.sub.0.sup.k=φ.sub.0.sup.k−1+α.sup.kp.sup.k
It can be shown that:
p.sup.k=−∇ε.sup.k+β.sup.kp.sup.k−1
where: ∇ε.sup.k is a gradient vector, of dimension (N.sub.pix, 1), each term of which represents a variation of the indicator ε.sup.k as a function of each of the degrees of freedom of the unknowns of the problem, that is, the terms of the vector φ.sub.0.sup.k; p.sup.k−1 is a direction vector established during the previous iteration; β.sup.k is a scale factor applied to the direction vector p.sup.k−1.
(64) Each term ∇ε.sup.k(x, y) of the gradient vector ∇ε is such that
(65)
where Im represents the imaginary operator and r′ represents a coordinate (x, y) in the detection plane.
(66) The scale factor β.sup.k may be expressed such that:
(67)
(68) The step α.sup.k may vary according to the iterations, for example, between 0.03 during the first iterations and 0.0005 during the last iterations.
(69) The updating equation makes it possible to obtain an adjustment of the vector φ.sub.0.sup.k, resulting in an iterative updating of the phase φ.sub.0.sup.k(x, y) in each pixel of the complex image A.sub.0.sup.k. (This complex image A.sub.0.sup.k, in the detection plane, is then updated by these new phase values associated with each pixel.
(70) Step 125: Reiteration or exiting from the algorithm.
(71) If no criterion of convergence has been achieved, step 125 consists in reiterating the algorithm, by a new iteration of steps 121 to 125 based on the complex image A.sub.0.sup.k updated during step 124. The criterion of convergence may be a predetermined number K of iterations, or a minimum value of the gradient ∇ε.sup.k of the indicator, or a difference considered to be negligible between two consecutive phase vectors φ.sub.0.sup.k−1, φ.sub.0.sup.k. When the criterion of convergence is achieved, one will have an estimation considered to be correct of a complex image of the sample in the detection plane P.sub.0.
(72) Step 126: Obtaining of the complex image in the reconstruction plane.
(73) At the end of the last iteration, the method involves a propagation of the complex image A.sub.0.sup.k resulting from the last iteration in the reconstruction plane P.sup.z, so as to obtain a complex image in the reconstruction plane A.sub.z=A.sub.z.sup.k.
(74) Experimental Trials
(75) During the course of a first series of trials, the method previously described was used with diluted blood samples in order to count the red blood cells.
(76) In this first series of trials, human blood was diluted in a sphering reagent, enabling a modification of the surface tension of the red blood cells so as to render them spherical. A preliminary calibration was done, with the aid of a reference machine, in order to produce, in a way familiar to the person skilled in the art, a recalibration between the concentration of the blood and the number of red blood cells counted. This recalibration takes into account the dilution factor, the thickness of the chamber, and the surface of the sample exposed to the image sensor.
(77) The protocol for preparation of the sample was as follows: dilution to 1/600 in a sphering reagent; taking a sample of 10 μl of diluted blood and injection into the fluidic chamber arranged opposite the image sensor;
(78) Each sample was the subject of a counting by an ABX Pentra DX 120 machine, used as the reference method.
(79) 80 samples were measured.
(80) A linear regression was established for each of these figures, the expression for each regression being respectively:
y=1.017x−0.04(r.sup.2=0.98);
y=1.019x−0.07(r.sup.2=0.98);
y=1.017x−0.07(r.sup.2=0.98);
y=1.027x−0.11(r.sup.2=0.98).
(81) The coefficient r.sup.2 is the coefficient of determination associated with each linear regression of Passing-Bablok type.
(82) During a second series of trials, we counted the white blood cells in samples of human blood, after lysis of the red blood cells by adding a lysis reagent. The image segmentation was done on the basis of a thresholding based on a maximum entropy criterion. In order to take account of the density of the white blood cells, the dilution factor used was equal to 10.
(83) The method was likewise tested during the course of a third series of trials, using glass beads of diameter 5 μm. These were reference beads of Bangs laboratories—SS06N, diluted to 1/2000 in a saline buffer of PBS type. The acquisition was done using a light source of laser type, as previously described. One observes the detection of a cluster containing two beads.
(84) In the case of red blood cells, the inventors believe it to be preferable for the surface density of particles of interest being counted to be less than 2000 particles per mm.sup.2 of surface of detection. Beyond this, the particles are too close together, which degrades the signal to noise ratio of each region of interest. Thus, it is advantageous to estimate, for the different cases in question, a maximum surface density and to then adapt the dilution being applied to the sample.
(85) One of the intended applications is the counting of blood cells, which may be red blood cells, white blood cells, or platelets. The inventors have discovered that in the case of white blood cells or red blood cells, during step 120, it is preferable for the reconstructed image I.sub.z to represent the modulus of the complex amplitude A.sub.z in the reconstruction plane P.sub.z. The definition of the edges of the regions of interest, corresponding to the cells, is more distinct. When the particles of interest are platelets, on the other hand, it seems preferable for the reconstructed image I.sub.z to represent the phase of the complex amplitude A.sub.z in the reconstruction plane P.sub.z.
(86) The inventors carried out experimental trials. After these trials, they defined the variants set forth below, each variant being able to be combined in whole or in part with the previously described steps.
(87) According to a first variant, one carries out the steps 100, 110 and 120 as previously described. During the course of step 120, the reconstruction distance taken into account corresponds to a focusing distance determined by carrying out a digital focusing algorithm, as previously described. Step 120 also involves a determination of a criterion of validity of the digital focusing. For this, the image acquired during step 110 is virtually partitioned into p zones, p being an integer which is strictly positive. Typically, p is between 1 and 100, such as 9. A digital focusing algorithm is applied in each zone so as to obtain, for each zone, a focusing distance d.sub.p. One then determines a dispersion indicator δ(d.sub.p) of the focusing distances.
(88) If the dispersion indicator δ(d.sub.p) exceeds a predetermined maximum dispersion indicator δ.sub.max, the digital focusing is invalidated. In this case, the focusing distance is considered to be equal to a mean focusing distance d.sub.mean established on the basis of a number k of previously analyzed samples, for which the focusing distance d.sub.k was memorized. For example, the mean focusing distance d.sub.mean may be calculated by a moving average, considering the last k samples analyzed, k being an integer, for example between 2 and 100. d.sub.mean may also be a median of the respective focusing distances d.sub.k of the last k samples analyzed.
(89) If the dispersion indicator δ(d.sub.p) is less than the maximum dispersion indicator δ.sub.max, the focusing distance is considered to be the mean
(90) The dispersion indicator δ(d.sub.p) may be a difference between a minimum focusing distance d.sub.p-min and a maximum focusing distance d.sub.p-max: δ(d.sub.p)=d.sub.p-max−d.sub.p-min It may also comprise a standard deviation or a variance of the focusing distances d.sub.p.
(91) Furthermore, in each zone partitioning the acquired image, the digital focusing algorithm is established by considering a lower focusing limit d.sub.inf and an upper focusing limit d.sub.sup. After the application of the digital focusing algorithm to each zone, if for at least one zone the focusing distance determined is less than or equal to the limit d.sub.inf, or greater than or equal to the limit d.sub.sup, the focusing distance is considered as being equal to the mean (or median) focusing distance d.sub.mean described above. In fact, if the focusing algorithm provides at least one focusing distance d.sub.p less than or equal to d.sub.inf or greater than or equal to d.sub.sup, it is probable that the focusing distance is in error. Thus, the focusing distance is replaced by the distance d.sub.mean, the latter being considered to be more probable.
(92) According to a second variant, in the course of step 130 the segmentation is done by an Otsu thresholding. The value of the segmentation threshold so obtained is compared to the mode of the histogram of the reconstructed image I.sub.z, whether a phase image or an image of the modulus, or a combination of these, of the reconstructed image. The mode of the histogram corresponds to a level of intensity of pixels corresponding to the maximum of the histogram of the reconstructed image. If the segmentation threshold determined by the Otsu thresholding is greater than 90% of the mode of the histogram, the segmentation threshold applied is arbitrarily set at a certain percentage of the mode of the histogram, for example, between 80% and 85% of the mode of the histogram. In fact, when the Otsu threshold is close to the mode of the histogram, the segmentation of the reconstructed image is generally not satisfactory, since zones of noise are considered as being one or more regions of interest. In these conditions, the fact of imposing a segmentation threshold between 80% and 85% of the mode of the histogram makes it possible to avoid this pitfall.
(93) According to a third variant, one determines a criterion of quality Q of the image I.sub.z* segmented during step 130. When the criterion of quality indicates a poor quality of the segmented image, the method is interrupted.
(94) The criterion of quality Q may correspond to an overall area of the regions of interest resulting from the segmentation performed during the course of step 130. Thus, one carries out steps 100, 110, 120 and 130 as previously described. During the course of step 130, one determines the total area A of the regions of interest determined by the segmentation, the latter being for example carried out by the application of Otsu thresholding. Then, if the total area A surpasses a predefined maximum area A.sub.max, the image acquired during step 110 is invalidated; if the total area A is less than the maximum area A.sub.max, but is situated close to the latter, one calculates an auxiliary criterion Q′, based on a signal to noise ratio SNR(I.sub.z*) of the segmented image I.sub.z*. By close to the maximum area A.sub.max is meant between the maximum area and a percentage of the maximum area A.sub.max, such as 80% or 90% of the maximum area.
(95) The maximum area A.sub.max may be determined by the person skilled in the art on the basis of experimental trials. The signal to noise ratio SNR(I.sub.z*) of the segmented image I.sub.z* corresponds to a level of intensity I.sub.ROI in the different regions of interest ROI segmented during step 130, normalized by an estimation of the noise N.sub.out outside of the regions of interest. The level of intensity I.sub.ROI may be estimated by a mean value μ or a median value med of the intensity in the different regions of interest ROI. The noise N.sub.out outside the regions of interest may be estimated by the standard deviation σ in the segmented image I.sub.z*, outside the regions of interest ROI. Thus:
(96)
(97) The signal to noise ratio so defined is then compared to a predetermined signal to noise ratio threshold SNR.sub.min. If the signal to noise ratio SNR(I.sub.z*) is greater than this threshold, the segmented image is considered to be valid. Otherwise, the segmented image is invalidated.
(98) According to such an embodiment, after step 130 one determines at least one criterion of quality Q with respect to the segmented image I.sub.z*. The criterion of quality may correspond to a comparison between an area A of all of the regions of interest segmented during the course of step 130 and a predefined maximum area A.sub.max. As a function of the comparison, the segmented image will be validated or invalidated. If the area A of all of the segmented regions of interest is less than the maximum area A.sub.max, close to the latter, an auxiliary criterion of quality Q′ will be determined. It may involve a comparison between the signal to noise ratio SNR(I.sub.z*) of the segmented image, as previously defined, and a minimum signal to noise ratio SNR.sub.min. As a function of this comparison, the segmented image will be validated or invalidated.
(99) This variant was applied to 305 samples, the particles being red blood cells. It was able to invalidate images corresponding to particular cases such as empty images, images in which the number of particles is too low or too high, or images resulting from a problem occurring during the acquisition.
(100) A fourth variant, described below, may relate more particularly to samples containing blood. During the course of experimental trials, the inventors considered that, as a general rule, when the particles of interest are red blood cells, it is advisable to provide, in step 140, a distinction between N types of regions of interest ROI.sub.n with 1≤n≤N. N is an integer generally between 3 and 5. It will be recalled that, as described in connection with step 140, the notation ROI.sub.n corresponds to a region of interest ROI containing n particles. The integer N corresponds to the maximum number of acceptable red blood cells in the same region of interest. When N=5, one distinguishes 5 types of regions of interest, corresponding to singlets (n=1), doublets (n=2), triplets (n=3), quadruplets (n=4) and quintuplets (n=5). It is considered that the regions of interest whose size exceeds an upper threshold S.sub.max, do not correspond to a cluster of red blood cells. Thus, step 140 involves taking into account a maximum size S.sub.max, and a rejection of a region of interest when its size is greater than the maximum size S.sub.max The maximum size S.sub.max is a size greater than the size corresponding to the region of interest ROI.sub.n=N. The maximum size S.sub.max is parametrizable, and depends on the experimental conditions, for example, the fluidic chamber used. The regions of interest whose size surpasses the maximum size S.sub.max are rejected. It has been found that, when the analyzed sample contains blood lacking in agglutinins, the number of regions of interest so rejected is less than 10. The number of regions of interest rejected can be determined by comparing: the number of regions of interest at the end of the segmentation step 130; the number of regions of interest at the end of the selection step 140.
(101) If the number of regions of interest rejected is less than or equal to a predefined level of acceptance L.sub.max, for example being less than or equal to 10, the count of the particles of interest is validated. The level of acceptance L.sub.max may be adjusted, by the person skilled in the art, during the course of an experimental learning phase. If the number of regions of interest rejected surpasses the level of acceptance L.sub.max, the count of the particles of interest is invalidated and this invalidation is reported. Such a situation may correspond, for example, to blood containing cold agglutinins. Owing to this, the blood contains many clusters of red blood cells, having a number of red blood cells greater than the maximum acceptable number of red blood cells N. In such a case, the inventors have observed that the number of regions of interest rejected might reach several tens, or even several hundreds. Given the number of regions of interest rejected, the count of the red blood cells is highly underestimated.
(102) During the course of experimental trials on six samples containing blood including cold agglutinins, and considering N=5, the inventors measured a number of rejected regions of interest respectively equal to 325, 379, 490, 508, 327 and 423. Given the number of rejected regions of interest, surpassing the level of acceptance L.sub.max, the count of the red blood cells was thus invalidated.
(103) According to such an embodiment, after step 130 one determines the number of regions of interest resulting from the segmentation of the image. Step 140 then involves: taking into account a maximum size S.sub.max of regions of interest, the regions of interest whose size exceeds the maximum size being rejected; a determination of the number of regions of interest rejected; taking into account a level of acceptance L.sub.max and a comparison between the number of regions of interest rejected and the level of acceptance, as a function of which the counting of the particles of interest is validated or invalidated.
(104) The level of acceptance L.sub.max may be defined by the person skilled in the art during the course of experimental trials.
(105) The invention may be used for the counting of particles of interest in the blood, but also in other bodily fluids, such as urine, saliva, sperm, etc. It may also be used in the quantification of microorganisms, such as bacteria or bacterial colonies. Beyond applications involving biology or healthcare, the invention may be used for inspection of samples in industrial fields or environmental monitoring, for example, in agriculture or the inspection of industrial fluids.