METHOD FOR OBSERVING A SAMPLE, BY CALCULATION OF A COMPLEX IMAGE

Abstract

The invention relates to a method for observing a sample (15), comprising the illumination of the sample using a light source (11) and the acquisition of an image (Io) of the sample using an image sensor (16), the sample being disposed between the image sensor and the light source. Iterative steps are applied to the acquired image (Io), also referred to as a hologram, comprising: a single iterative numerical propagation (h), such as to estimate a complex image (A) of the sample in a reconstruction plane (P10) or in a detection plane (P0), in which the image sensor extends. The complex image can be used for the characterisation of the sample.

Claims

1. A method for observing a sample comprising: a) illuminating the sample using a light source, configured to emit a light wave that propagates towards the sample; b) acquiring, using an image sensor, an image of the sample, the image being formed in a detection plane, the sample being placed between the light source and the image sensor, the acquired image being representative of a light wave transmitted by the sample under the effect of the illumination, the acquired image comprising interference patterns formed by interference between the light wave emitted by the light source and a light wave diffracted by the sample; and c) calculating a complex image representative of the sample; wherein calculating the complex image representative of the sample comprises: i) defining an initial image of the sample in the detection plane, on the basis of the image acquired by the image sensor; ii) determining a complex image of the sample in a reconstruction plane by applying a propagation operator to an initial image of the sample, which is the image defined in i), or to the image of the sample, in the detection plane, resulting from the preceding iteration; iii) calculating a noise indicator on the basis of the complex image obtained in ii), the noise indicator depending on a reconstruction noise affecting the complex image, the noise indicator being representative of fluctuations appearing in the complex image; iv) updating the image of the sample in the detection plane by adjusting phase values of the pixels of the d image, the adjustment being carried out depending on a variation in the noise indicator calculated in iii) with respect to the phase values; v) reiterating ii) to iv) until a convergence criterion is reached, so as to obtain a complex image of the sample in the detection plane, or in the reconstruction plane.

2. The method of claim 1, wherein iii) comprises: for various pixels, calculating a quantity associated with each pixel, depending on the value of the complex image determined in ii) for the pixel, or on a dimensional derivative of the complex image for the pixel; and combining the quantities calculated for different pixels, so as to obtain the noise indicator.

3. The method of claim 2, wherein the noise indicator is a norm of order lower than or equal to 1 calculated on the basis of the quantities respectively associated with each pixel.

4. The method of claim 2, wherein, in iii), the quantity associated with each pixel is calculated on the basis of a modulus of a dimensional derivative, for the pixel, of the complex image determined in substep ii).

5. The method of claim 3, wherein, in iii), the noise indicator is calculated on the basis of the value of the complex image in a plurality of pixels.

6. The method of claim 5, wherein: in i), the initial image of the sample is defined by normalizing the image acquired by the image sensor by an image representative of the light wave emitted by the light source; in iii), the quantity associated with each pixel is calculated depending on the value of the complex image determined in ii), for the pixel, subtracted from a strictly positive number.

7. The method of claim 2, wherein, in iii), the noise indicator is a sum, which is optionally weighted, of the quantity associated with each pixel of the complex image in the reconstruction plane.

8. The method of claim 1, wherein, in iv), the value of the phase of each pixel is adjusted by forming a phase vector, each term of which corresponding to the value of the phase of a pixel of the image of the sample in the detection plane, the phase vector being updated, in each iteration, so as to minimise, or to maximise, the noise indicator calculated in iii), the update being based on a gradient of the noise indicator, with respect to each term of the phase vector.

9. The method of claim 8, wherein the phase vector is updated by a minimization or maximization algorithm based on a gradient of the noise indicator with respect to each term of the phase vector.

10. The method of claim 1, wherein the reconstruction plane lies parallel to the detection plane.

11. The method of claim 1, wherein the reconstruction plane is a plane in which the sample lies.

12. The method of claim 1, wherein no magnifying optic is interposed between the sample and the image sensor.

13. The method of claim 1, further comprising d) characterizing the sample depending on the calculated complex image.

14. A device for observing a sample, the device comprising: a light source configured to emit an incident light wave that propagates towards the sample; an image sensor; a holder, configured to hold the sample between the light source and the image sensor; and a processor, configured to receive an image of the sample acquired by the image sensor and to implement the method of claim 1.

Description

FIGURES

[0041] FIG. 1 shows an example of a device according to the invention.

[0042] FIG. 2A illustrates the main steps of a method allowing a complex image of the sample in a reconstruction plane to be calculated.

[0043] FIGS. 2B, 2C, 2D and 2E respectively show:

[0044] a hologram acquired by the image sensor;

[0045] an image reconstructed in a reconstruction plane in a first iteration of the method shown in FIG. 2A;

[0046] an image showing a quantity associated with each pixel of the image shown in FIG. 20, and

[0047] a representation of a reconstructed image after a plurality of iterations of the method shown in FIG. 2A.

[0048] FIGS. 3A and 3B schematically show a profile of the modulus or of the phase of a complex image obtained by holographic reconstruction, in the presence of and without reconstruction noise, respectively.

[0049] FIG. 4 summarizes the operation of a method implementing the invention.

[0050] FIG. 5A is a hologram acquired by an image sensor, the sample comprising cells dispersed in a continuous aqueous phase. FIGS. 5B and 5C show the modulus and the phase of a complex image in the plane of the sample, this complex image being formed by implementing the invention, respectively.

[0051] FIG. 6A is a hologram acquired by an image sensor, the sample containing red blood cells dispersed in a continuous aqueous phase; and FIGS. 6B and 6C show the modulus and the phase of a complex image in the plane of the sample, this complex image being formed by implementing the invention, respectively.

DESCRIPTION OF PARTICULAR EMBODIMENTS

[0052] FIG. 1 shows an example of a device according to the invention. A light source 11 is able to emit a light wave 12, called the incident light wave, that propagates in the direction of a sample 10, along a propagation axis Z. The light wave is emitted in a spectral band including a wavelength . This wavelength may be a central wavelength of said spectral band.

[0053] The sample 10 is a sample that it is desired to characterize. It may in particular be a question of a medium 10a containing particles 10b. The particles 10b may be blood particles, for example red blood cells. It may also be a question of cells, microorganisms, for example bacteria or yeasts, microalgae, microbeads, or droplets that are insoluble in the liquid medium, for example lipid nanoparticles. Preferably, the particles 10b have a diameter, or are inscribed in a diameter, smaller than 1 mm, and preferably smaller than 100 m. It is a question of microparticles (diameter smaller than 1 mm) or of nanoparticles (diameter smaller than one micron). The medium 10a, in which the particles are immersed, may be a liquid medium, for example a liquid phase of a bodily liquid, a culture medium or a liquid sampled from the environment or from an industrial process. It may also be a question of a solid medium or a medium having the consistency of a gel, for example an agar substrate, propitious to the growth of bacterial colonies.

[0054] The sample 10 is contained in a fluidic chamber 15. The fluidic chamber 15 is for example a microcuvette, commonly used in point-of-care type devices, into which fluidic chamber 15 the sample 20 penetrates, for example via capillary action. The thickness e of the sample 10, along the propagation axis, typically varies between 20 m and 1 cm, and is preferably comprised between 50 m and 500 m, and is for example 150 m.

[0055] The sample lies in a plane P.sub.10, called the sample plane, that is perpendicular to the propagation axis. It is held on a holder 10s.

[0056] The distance D between the light source 11 and the sample 10 is preferably larger than 1 cm. It is preferably comprised between 2 and 30 cm. Preferably, the light source, seen by the sample, may be considered to be point-like. This means that its diameter (or its diagonal) is preferably smaller than one tenth, and better still one hundredth, of the distance between the sample and the light source. Thus, preferably, the light reaches the sample in the form of plane waves, or waves that may be considered as such.

[0057] The light source 11 may be a light-emitting diode or a laser diode. It may be associated with a diaphragm 18, or spatial filter. The aperture of the diaphragm is typically comprised between 5 m and 1 mm, and preferably between 50 m and 500 m. In this example, the diaphragm is supplied by Thorlabs under the reference P150S and its diameter is 150 m. The diaphragm may be replaced by an optical fiber, a first end of which is placed facing the light source 11 and a second end of which is placed facing the sample 10.

[0058] The device preferably includes a diffuser 17, placed between the light source 11 and the diaphragm 18. Use of such a diffuser allows constraints on the centering of the light source 11 with respect to the aperture of the diaphragm 18 to be relaxed. The function of such a diffuser is to spread the light beam, produced by an elementary light source 11.sub.i, (1i3) over a cone of angle , a being equal to 30 in the present case. Preferably, the diffusion angle varies between 10 and 80.

[0059] Preferably, the spectral emission band of the incident light wave 12 has a width smaller than 100 nm. By spectral bandwidth, what is meant is a full width at half maximum of said spectral band.

[0060] The sample 10 is placed between the light source 11 and an image sensor 16. The latter preferably lies parallel, or substantially parallel, to the plane in which the sample lies. The expression substantially parallel means that the two elements may not be rigorously parallel, an angular tolerance of a few degrees, smaller than 20 or 10 being acceptable.

[0061] The image sensor 16 is able to form an image, in a detection plane P.sub.0. In the example shown, it is a question of an image sensor comprising a matrix array of pixels, of CCD or CMOS type. CMOS sensors are the preferred type of sensor because the pixels are smaller in size, this allowing images the spatial resolution of which is more favorable to be acquired. The detection plane P.sub.0 preferably lies perpendicular to the propagation axis Z of the incident light wave 12.

[0062] The distance d between the sample 10 and the matrix array of pixels of the image sensor 16 is preferably comprised between 50 m and 2 cm, and preferably comprised between 100 m and 2 mm.

[0063] The absence of magnifying optics between the image sensor 16 and the sample 10 will be noted. This does not prevent focusing micro-lenses optionally being present near each pixel of the image sensor 16, the latter not having the function of magnifying the image acquired by the image sensor.

[0064] Under the effect of the incident light wave 12, the sample 10 may generate a diffracted wave, that is liable to produce, in the detection plane P.sub.0, interference, in particular with a portion of the incident light wave 12 that is transmitted by the sample. Moreover, the sample may absorb a portion of the incident light wave 12. Thus, the light wave 22, transmitted by the sample, and to which the image sensor 20 is exposed, may comprise: [0065] a component resulting from the diffraction of the incident light wave 12 by the sample; [0066] a component resulting from the absorption of the incident light wave 12 by the sample.

[0067] The light wave 22 could also be called the exposure light wave. A processor 20, for example a microprocessor, is able to process each image 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 to carry out the image-processing operations and calculations described in this description is stored. The processor may be coupled to a screen 24 allowing the display of images acquired by the image sensor 16 or calculated by the processor 20.

[0068] An image acquired on the image sensor 16, which image is also called a hologram, does not allow a sufficiently precise representation of the observed sample to be obtained. As described with reference to the prior art, it is possible to apply, to each image acquired by the image sensor, a propagation operator h, so as to calculate a quantity representative of the light wave 22 transmitted by the sample 10, and to which the image sensor 16 is exposed. Such a method, which is referred to as a holographic reconstruction method, in particular allows an image of the modulus or of the phase of this light wave 22, in a reconstruction plane parallel to the detection plane P.sub.0, and in particular in the plane P.sub.10 in which the sample lies, to be reconstructed. To do this, the image I.sub.0 acquired by the image sensor 16 and a propagation operator h are convoluted. It is then possible to reconstruct a complex expression A for the light wave 22 at every point of spatial coordinates (x,y,z), and in particular in a reconstruction plane P.sub.z located at a distance |z| from the image sensor 16, this reconstruction plane possibly being the plane P.sub.10 of the sample. The complex expression A is a complex quantity the argument and the modulus of which are representative of the phase and of the intensity of the light wave 22 to which the image sensor 16 is exposed, respectively. Convoluting the image I.sub.0 and the propagation operator h allows a complex image A.sub.z representing a spatial distribution of the complex expression A in a plane P.sub.z, called the reconstruction plane lying at a coordinate z from the detection plane P.sub.0, to be obtained. In this example, the equation of the detection plane P.sub.0 is z=0. This complex image corresponds to a complex image of the sample in the reconstruction plane P.sub.z. It also represents a two-dimensional spatial distribution of the optical properties of the wave 22 to which the image sensor 16 is exposed.

[0069] The function of the propagation operator h is to describe the propagation of light between the image sensor 16 and a point of coordinates (x,y,z), located at a distance |z| from the latter. It is then possible to determine the modulus M(x,y,z) and/or the phase (x,y,z) of the light wave 22, at this distance |z|, which is called the reconstruction distance, with:

M(x,y,z)=abs[A(x,y,z)](1)

(x,y,z)=arg[A(x,y,z)](2)

[0070] The operators abs and arg are the modulus and argument, respectively.

[0071] In other words, the complex amplitude A of the light wave 22 at every point of spatial coordinates (x,y,z) is such that: A (x,y,z)=M(x,y,z)e.sup.j(x,y,z) (3).

[0072] However, such a reconstruction is accompanied by reconstruction noise that may be substantial, because the propagation is carried out on the basis of an image I.sub.0 containing no phase information.

[0073] An algorithm allowing a complex image, denoted A.sub.z, representing a spatial distribution of the complex amplitude A(x,y,z) at a distance |z|, to be obtained will now be described with reference to FIG. 2A, the images obtained in certain steps being illustrated in FIGS. 2B to 2E. Such an image is representative of the light wave 22, to which the image sensor 16 is exposed, and in particular allows the modulus or the phase of this light wave to be accessed.

Step 100: Image Acquisition:

[0074] In this step, the image sensor 16 acquires an image I.sub.0 of the sample 16, and more precisely of the light wave 22 transmitted by the latter, to which the image sensor is exposed. Such an image, or hologram, is shown in FIG. 2B.

[0075] This image was produced using a sample 10 containing Chinese hamster ovary (CHO) cells submerged in a saline buffer, the sample being contained in a fluidic chamber of 100 m thickness placed as a distance d of 1500 m from a CMOS sensor. The sample was illuminated with a light-emitting diode 11 the spectral emission band of which is centered on a wavelength of 450 nm and located at a distance D=8 cm from the sample.

Step 110: Initialization:

[0076] In this step, an initial image A.sub.0.sup.k=0 of the sample is defined on the basis of the image I.sub.0 acquired by the image sensor 16. This step is an initialisation of the iterative algorithm described below with reference to the steps 120 to 180, the exponent k giving the rank of each iteration. 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)}. It may also be obtained by normalising the image I.sub.0 with a term representative of the intensity of the light wave 12 incident on the sample 16. The latter may be: [0077] the square root of an average I.sub.0 of the image I.sub.0, in which case each pixel I.sub.0 (x,y) of the acquired image is divided by said average, such that

[00001]$\begin{array}{cc}{M}_0^{k=0}=\sqrt{\frac{I_0\left(x,y\right)}{\stackrel{\_}{I_0}}}& \left(4\right)\end{array}$ [0078] an image I.sub.12 without sample, acquired by the image sensor 16 in the absence of sample between the light source 11 and the image sensor, in which case the value I.sub.0(x,y) of each pixel of the acquired image of the sample is divided by the value I.sub.12(x,y) of each pixel of the image without sample:

[00002]$M_0^{k=0}=\sqrt{\frac{I_0\left(x,y\right)}{I_{12}\left(x,y\right)}}$

(4) [0079] an average I.sub.12 of said image without sample, in which case each pixel I.sub.0(x,y) of the acquired image is divided by said average:

[00003]$\begin{array}{cc}{M}_0^{k=0}=\sqrt{\frac{I_0\left(x,y\right)}{\stackrel{\_}{I_{12}}}}& \left(4^{}\right)\end{array}$

[0080] The phase .sub.0.sup.k=0 of the initial image A.sub.0.sup.k=0 is either considered to be zero for each pixel (x,y), or preset to an arbitrary value. Specifically, the initial image Ar.sub.0.sup.k=0 results directly from the image I.sub.0 acquired by the image sensor 16. However, the latter image contains no information on the phase of the light wave 22 transmitted by the sample 10 the image sensor 16 being sensitive only to the intensity of this light wave.

Step 120: Propagation:

[0081] In this step, the image A.sub.0.sup.k1 obtained in the plane of the sample is propagated to a reconstruction plane P.sub.z, by applying a propagation operator such as described above, so as to obtain a complex image A.sub.z.sup.k, which is representative of the sample, in the reconstruction plane P.sub.z. The wording complex image expresses the fact that each term of this image is a complex quantity. The propagation is achieved by convoluting the image A.sub.z.sup.k1 with the propagation operator h.sub.z, such that:

A.sub.k.sup.z=A.sub.0.sup.k1*h.sub.z,(5),

the symbol * representing a convolution. The index z represents the fact that the propagation is carried out in a direction opposite to the propagation axis Z. Back propagation is spoken of.

[0082] The propagation operator is for example the Fresnel-Helmholtz function, such that:

[00004]$\begin{array}{cc}h\left(x,y,z\right)=\frac{1}{j.Math..Math..Math..Math.z}.Math.e^{j.Math..Math.2.Math..Math..Math..Math.\frac{z}{}}.Math.\exp \left(j.Math..Math..Math..Math.\frac{x^2+y^2}{.Math..Math.z}\right)& \left(6\right)\end{array}$

[0083] The convolution is generally carried out in the frequency domain, in which it equates to a product, in which case the Fourier transform of this operator is used, the latter being:

[00005]$\begin{array}{cc}H\left(,v,z\right)=e^{j.Math..Math.2.Math..Math..Math.\frac{z}{}}.Math.\exp \left(-.Math..Math.z\left({}^2+v^2\right)\right)& \left(6^{}\right)\end{array}$

where is the central wavelength of the spectral emission band of the light source 11.

[0084] Thus,

[00006]$\begin{array}{cc}{A}_z^k\left(x,y\right)=A_z^k\left(r\right)=-\frac{1}{j.Math..Math..Math..Math.z}.Math.e^{-j.Math..Math.2.Math..Math..Math..Math.\frac{z}{}}.Math.A_0^{k-1}\left(r^{}\right).Math.\exp -\left(j.Math..Math..Math..Math.\frac{{\left(r-r^{}\right)}^2}{.Math..Math.z}\right).Math.{\mathrm{dr}}^{}& \left(7\right)\end{array}$

[0085] where r and r are radial coordinates, i.e. in the reconstruction plane P.sub.z and in the detection plane P.sub.0, respectively.

[0086] In the first iteration (k=1), A.sub.0.sup.k=0 is the initial image determined in step 110. In following iterations, A.sub.0.sup.k1, is the complex image updated in the preceding iteration, in the detection plane P.sub.0.

[0087] The reconstruction plane P.sub.z is a plane that is distant from the detection plane P.sub.0, and preferably parallel to the latter. Preferably, the reconstruction plane P.sub.z is a plane P.sub.10 in which the sample 10 lies. Specifically, an image reconstructed in this plane allows a generally high spatial resolution to be obtained. It may also be a question of another plane, located at a nonzero distance from the detection plane, and preferably parallel to the latter, for example a plane lying between the image sensor 16 and the sample 10.

[0088] FIG. 2C shows the modulus of an image A.sub.z.sup.k=1, reconstructed at a distance of 1440 m from the detection plane P.sub.0 by applying the propagation operator defined above to the hologram of FIG. 2B. This image shows the complex image, in the reconstruction plane, established in the first iteration.

Step 130: Calculating a Quantity for a Plurality of Pixels of the Complex Image:

[0089] In this step, 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 with each of these pixels, is calculated. This quantity depends on the value A.sub.z.sup.k(x,y) of the image A.sub.z.sup.k, or on its modulus, in the pixel (x,y) for which it is calculated. It may also depend on a dimensional derivative of the image in this pixel, for example the modulus of a dimensional derivative of this image.

[0090] In this example, the quantity .sup.k(x,y) associated with each pixel (x,y) is based on the modulus of a dimensional derivative, such that:

[00007]$\begin{array}{cc}{\mathrm{.Math.}}^k\left(x,y\right)=\sqrt{{.Math.\frac{A_z^k\left(x,y\right)}{x}.Math.}^2+{.Math.\frac{A_z^k\left(x,y\right)}{y}.Math.}^2}& \left(8\right)\end{array}$

[0091] Since the image is discretized into pixels, the derivative operators may be replaced by Sobel operators, such that:

[00008]$\begin{array}{cc}{\mathrm{.Math.}}^k\left(x,y\right)=\sqrt{\left(S_x*A_z^k\left(x,y\right)\right).Math.{\left(S_x*A_z^k\left(x,y\right)\right)}^{*}+\left(S_y*A_z^k\left(x,y\right)\right).Math.{\left(S_y*A_z^k\left(x,y\right)\right)}^{*}}& \left(9\right)\end{array}$

where: [0092] ( )* represents the complex-conjugate operator; and [0093] S.sub.x and S.sub.y are Sobel operators along two orthogonal axes X and Y of the reconstruction plane P.sub.z.

[0094] In this example,

[00009]$\begin{array}{cc}{S}_x=\left[\begin{array}{ccc}1& 0& -1\\ 2& 0& -2\\ 1& 0& -1\end{array}\right]& \left(10\right)\end{array}$

and S.sub.y is the transpose matrix of S.sub.x.

[0095] FIG. 2D shows, in the form of an image, the value of the quantity .sup.k(x,y) in each pixel of the image A.sub.z.sup.k=1 shown in FIG. 2C.

Step 140: establishing a noise indicator associated with the image A.sub.z.sup.k:

[0096] In step 130, quantities .sup.k(x,y) are calculated for a plurality of pixels of the complex image A.sub.z.sup.k. These quantities may form a vector E.sup.k, the terms of which are the quantities .sup.k(x,y) associated with each pixel (x,y).

[0097] In this step, an indicator, called the noise indicator, is calculated on the basis of a norm of the vector E.sup.k. Generally, an order is associated with a norm, such that the norm x.sub.p of order p of a vector x of dimension n of coordinates (x.sub.1, x.sub.2, . . . x.sub.n,) is such that:

x.sub.p=(.sub.i=1.sup.n|x.sub.i|x.sub.i|.sup.p).sup.1/p, with p0.(12)

[0098] In the present case, a norm of order 1, is used, in other words p=1. Specifically, the inventors have come to the conclusion that use of a norm of order 1 or of order lower than or equal to 1, is particularly suited to this algorithm, for the reasons explained below with reference to FIGS. 3A and 3B.

[0099] In this step, the quantity .sup.k (x,y) calculated on the basis of the complex image A.sub.z.sup.k for each pixel (x,y) of the latter, is summed so as to form a noise indicator .sup.k associated with the complex image A.sub.z.sup.k.

Thus, .sup.k=.sub.(x,y).sup.k(x,y)(15).

[0100] This noise indicator .sup.k corresponds to a total variation norm of the complex image A.sub.z.sup.k.

[0101] With reference to the image in FIG. 2D, the noise indicator .sup.k=1 is obtained, in the first iteration, by summing the value of the pixels of this image.

[0102] Alternatively to a norm of order 1, a weighted sum of the quantities .sup.k (x,y), or another arithmetic combination, is also envisionable.

[0103] Because of the use of a norm of order 1, or of order lower than or equal to 1, the value of the noise indicator .sup.k decreases as the complex image A.sub.z.sup.k becomes more representative of the sample. Specifically, in the first iterations, the value of the phase .sub.0.sup.k(x,y), in the each pixel (x,y) of the image A.sub.0.sup.k, is poorly estimated. The propagation of the image of the sample from the detection plane P.sub.0 to the reconstruction plane P.sub.z is thus accompanied by a substantial amount of reconstruction noise, as mentioned with reference to the prior art. This reconstruction noise takes the form of fluctuations that appear in the reconstructed image. Because of these fluctuations, a noise indicator .sup.k, such as defined above, increases as the contribution of the reconstruction noise, to the reconstructed image, increases. Specifically, fluctuations due to the reconstruction noise tend to increase the value of this indicator.

[0104] FIGS. 3A and 3B schematically show a radial profile of the modulus (or of a phase) of a reconstructed image, that is affected by a large amount and a small amount of reconstruction noise, respectively. A sample containing a dispersion of particles 10b in a transparent homogenous medium 10a is considered here. The schematic profiles contain two large fluctuations, each of these fluctuations being representative of a particle 10b. The profile of FIG. 3A also includes high-frequency fluctuations of lower amplitude that are representative of reconstruction noise. The noise indicator .sup.k, such as defined above, is higher for FIG. 3A than for FIG. 3B. An indicator .sup.k based on a norm of order higher than 1 could also be used, but such a norm tends to attenuate low-amplitude fluctuations, which are representative of reconstruction noise, with respect to large fluctuations, which are representative of the sample. In contrast, a norm of order 1, or of order lower than 1, does not attenuate small fluctuations with respect to large fluctuations. This is why the inventors prefer an indicator of reconstruction noise .sup.k based on a norm of order 1 or lower than 1.

[0105] An important aspect of this step consists in determining, in the detection plane P.sub.0, for each pixel of the image A.sub.0.sup.k of the sample, phase values .sub.0.sup.k(x,y) allowing, in a following iteration, a reconstructed image A.sub.z.sup.k+1 to be obtained the indicator .sup.k+1 of which is lower than the indicator .sup.k.

[0106] In the first iteration, as explained above, relevant information is only available on the intensity of the light wave 22 and not on its phase. The first image A.sub.z.sup.k=1 reconstructed in the reconstruction plane P.sub.z is therefore affected by a substantial amount of reconstruction noise, because of the absence of relevant information as to the phase of the light wave 22 in the detection plane P.sub.0. Therefore, the indicator .sup.k=1 is high. In following iterations, the algorithm gradually adjusts the phase .sub.0.sup.k(x,y) in the detection plane P.sub.0, so as to gradually minimise the indicator .sup.k.

[0107] The image A.sub.0.sup.k in the detection plane is representative of the light wave 22 in the detection plane P.sub.0, both from the point of view of its intensity and of its phase. The steps 120 to 160 aim to establish, iteratively, for each pixel of the image A.sub.0.sup.k, the value of the phase .sub.0.sup.k(x,y) that minimizes the indicator .sup.k, the latter being obtained A.sub.z.sup.k from the image obtained by propagation of the image A.sub.0.sup.k1 to the reconstruction plane P.sub.z.

[0108] The minimization algorithm may be a gradient descent algorithm or a conjugated gradient descent algorithm, the latter being described below.

Step 150: Adjusting the Value of the Phase in the Detection Plane:

[0109] Step 150 aims to determine a value for the phase .sub.0.sup.k(x,y) of each pixel of the complex image A.sub.0.sup.k so as to minimise the indicator .sup.k+1 resulting from a propagation of the complex image A.sub.0.sup.k to the reconstruction plane P.sub.z, in the following iteration k+1.

[0110] To do 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 is the number of pixels in question. This vector is updated in each iteration, using the following updating expression:

.sub.0.sup.k(x,y)=.sub.0.sup.k1(x,y)+.sup.kp.sup.k(x,y)(16) where:

[0111] .sup.k is a scalar, referred to as the step size, and representing a distance; [0112] 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.

[0113] This equation may be expressed in vector form, as follows:

.sub.0.sup.k=.sub.0.sup.k1+.sup.kp.sup.k(16)

[0114] It is possible to show that:

p.sup.k=.sup.k+.sup.kp.sup.k1(17)

where: [0115] .sup.k is a gradient vector, of dimension (N.sub.pix, 1), each term of which represents a variation in the indicator .sup.k as a function of each of the degrees of freedom, forming the unknowns of the problem, i.e. the terms of the vector .sub.0.sup.k; [0116] p.sup.k1 is a direction vector established in the preceding phase; [0117] .sup.k is a scale factor applied to the direction vector p.sup.k1.

[0118] Each term .sup.k(x,y) of the gradient vector .sup.k, is such that

[00010]$\begin{array}{cc}{\mathrm{.Math.}}^k\left(r^{}\right)=\frac{{\mathrm{.Math.}}^k}{{}_0^k\left(r^{}\right)}=-\mathrm{Im}\left(A_0^{k^{*}}\left(r^{}\right).Math.\left(\left(S_x*\frac{S_x*A_z^k}{{\mathrm{.Math.}}^k}+\frac{S_y*A_z^k}{{\mathrm{.Math.}}^k}.\right)*h_z\right).Math.\left(r^{}\right)\right)& \left(18\right)\end{array}$

where Im is the imaginary-part operator and r is a coordinate (x,y) in the detection plane.

[0119] The scale factor .sup.k is a scalar, it may be expressed such that:

[00011]$\begin{array}{cc}{}^{\left(k\right)}=\frac{{\mathrm{.Math.}}^{\left(k\right)}.Math.{\mathrm{.Math.}}^{\left(k\right)}}{{\mathrm{.Math.}}^{\left(k-1\right)}.Math.{\mathrm{.Math.}}^{\left(k-1\right)}},& \left(19\right)\end{array}$

the symbol .Math. representing a scalar product.

[0120] The step size .sup.k may vary depending on the iteration, for example between 0.03 in the first iterations and 0.0005 in the last iterations.

[0121] The updating equation allows an adjustment of the vector .sub.0.sup.k to be obtained, thereby leading to an iterative update of the phase .sub.0.sup.k(x,y) of 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 with these new values of the phase associated with each pixel. It will be noted that the modulus of the complex image A.sub.0.sup.k is not modified, the latter being determined on the basis of the image acquired by the image sensor, such that M.sub.0.sup.k(x,y)=M.sub.0.sup.k=0(x,y).

Step 160: Reiteration or Termination of the Algorithm:

[0122] Provided that a convergence criterion has not been met, step 160 consists in reiterating the algorithm, with a new iteration of steps 120 to 160, on the basis of the complex image A.sub.0.sup.k updated in step 150.

[0123] The convergent criterion may be a preset 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.k1, .sub.0.sup.k. When the convergence criterion has been met, an estimation that is considered to be correct of a complex image of the sample, in the detection plane P.sub.0 or in the reconstruction plane 4, is considered to have been obtained.

[0124] The method may comprise an ultimate step 170 of reconstructing a complex image A.sub.z of the sample in the reconstruction plane P.sub.z, by applying a propagation operator to the image A.sub.0.sup.k of the sample in the detection plane obtained in the last iteration.

[0125] FIG. 2E shows an image of the modulus M.sub.Z.sup.k30 of each pixel of the reference complex image A.sub.Z.sup.k=30 obtained in a reconstruction plane P.sub.z at the end of 30 iterations. This image may be compared to FIG. 2C, showing a similar image A.sub.Z.sup.k=1 obtained in the first iteration. A clear decrease in reconstruction noise, in particular between each pixel, may be seen.

[0126] In the example that has been given, the quantity .sup.k(x,y), associated with each pixel, and implemented in step 130, is based on a dimensional derivative calculated for each pixel (x,y) of the image A.sub.Z.sup.k. According to one variant, the initial image A.sub.0.sup.k=0 is normalised, as described above, by a scalar or an image representative of the incident wave 12. In this way, for each pixel, the modulus of the image of the sample, in the detection plane or in the reconstruction plane, is lower than or equal to 1. The quantity .sup.k(x,y) associated with each pixel, in step 130, is a modulus of a difference between the image A.sub.z.sup.k, in each pixel, and the value 1. Such a quantity may be obtained using 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|(20)

and, in step 140, .sup.k=.sub.(x,y).sup.k(x,y) (21), this corresponding, in a non-discretized form, to .sup.k(A.sub.z.sup.k)=|A.sub.z.sup.k1|=d{right arrow over (r)}.Math.{square root over ((A.sub.z.sup.k(r)1)(A.sub.z.sup.k(r)1)*)} (22), r being a radial coordinate in the reconstruction plane.

[0127] The noise indicator is yet again a norm of order 1 of a vector .sup.k each term of which is the modulus .sup.k(x,y) calculated for each pixel.

[0128] It is possible to show that the gradient of this noise indicator .sup.k, with respect to the phase vector, is such that:

[00012]$\begin{array}{cc}s.Math..Math.{\mathrm{.Math.}}^k\left(r^{}\right)=\frac{{\mathrm{.Math.}}^k}{{}_0^k\left(r^{}\right)}=\mathrm{Im}\left(A_0^{k^{*}}\right).Math.\left(r^{}\right).Math.\left(\frac{\left(A_z^k-1\right)}{.Math.A_z^k-1.Math.}*h_z\right).Math.\left(r^{}\right)& \left(23\right)\end{array}$

r being a radial coordinate in the detection plane.

[0129] Just like the total variation norm described above with reference to expression (15), the use of such an indicator is suited to a sample containing particles 10b dispersed in a homogenous medium 10a. During the gradient descent algorithm, this indicator tends to decrease the number of pixels the modulus of which is not equal to 1 in zones that are distributed discretely in the image of the sample, these zones corresponding to the particles 10b of the sample.

[0130] According to another variant, in step 130, the norm is such that:

.sup.k(x,y)=((A.sub.z.sup.k(x,y))(A.sub.z.sup.k(x,y))*1)=|A.sub.z.sup.k(x,y)|.sup.21(25)

Then, in step 140, .sup.k=.sub.(x,y)(.sup.k(x,y)).sup.2 (26), which corresponds, in a non-discretized form, to .sup.k=(|A.sub.z.sup.k|.sup.21).sup.2 (27)

[0131] Just as in the preceding embodiments, step 130 includes calculating a quantity .sup.k(x,y) associated with each pixel, based on a modulus of the complex image A.sub.z.sup.k, then calculating a noise indicator associated with the complex image A.sub.z.sup.k based on a norm. According to this variant, it is a question of a norm of order 2.

[0132] It is possible to show that the gradient of this indicator, with respect to the phase vector, is such that:

[00013]$\begin{array}{cc}{\mathrm{.Math.}}^k\left(r^{}\right)=\frac{{\mathrm{.Math.}}^k}{{}_0^k\left(r^{}\right)}=2.Math..Math.\mathrm{Im}.Math..Math.A_0^{k^{*}}\left(\left({.Math.A_z^k.Math.}^2-1\right).Math.A_z^k\right)*h_z\left(r^{}\right)& \left(28\right)\end{array}$

[0133] According to another variant, in step 130, the quantity associated with each pixel is such that

.sup.k(x,y)=.sup.k(x,y)=((A.sub.z.sup.k(x,y))(A.sub.z.sup.k(x,y))*1)=|A.sub.z.sup.k(x,y)|.sup.21(30)

[0134] Then, in step 140, .sup.k=.sub.(x,y).sup.k(x,y) (31), which corresponds, in a non-discretized form, to .sup.k=A.sub.z.sup.k|.sup.21| (32). The quantity associated with each pixel is identical to the preceding variant (see equation (25)), but the noise indicator associated with the image is calculated using a norm of order 1.

[0135] According to another variant, in step 130, the quantity associated with each pixel is such that

.sup.k(x,y)=I{square root over ((A.sub.z.sup.k(x,y))(A.sub.z.sup.k(x,y))*)}1|=|A.sub.z.sup.k(x,y)|1(35)

[0136] Then, in step 140, .sup.k=.sub.(x,y).sup.k(x,y) (36), which corresponds, in a non-non-discretized form, to .sup.k=A.sub.z.sup.k|1| (37)

[0137] Thus, whatever the embodiment, the noise indicator .sup.k associated with a complex image A.sub.z.sup.k, may be obtained by: [0138] calculating a quantity for a plurality of pixels of the image, based on a value of the latter, a modulus of the latter or a dimensional derivative of the latter; and [0139] combining said quantities in the form of a norm, and preferably of a norm of order lower than 1.

[0140] FIG. 4 summarizes the main steps of the algorithm described above: on the basis of an image I.sub.0 acquired by the image sensor 16, an initial image A.sub.0.sup.k=0 is formed. In each iteration k, a complex image A.sub.z.sup.k, representing the light wave 22 in a reconstruction plane P.sub.z, is established by numerical propagation of an image A.sub.0.sup.k representing the light wave 22 in the detection plane P.sub.0. A noise indicator .sup.k is associated with the image A.sub.z.sup.k. Its gradient .sup.k, as a function of the phase .sub.0.sup.k of the light wave 22, in the detection plane P.sub.0, is calculated, on the basis of which said phase of the light wave 22 in the detection plane, is updated. This update allows a new complex image A.sub.0.sup.k+1 to be formed in the detection plane P.sub.0, on the basis of which image a new iteration may be carried out. A notable aspect of the invention is that the phase .sub.0.sup.k is updated in the detection plane P.sub.0 using an optimisation algorithm based on the gradient of the indicator as a function of said phase.

[0141] In the embodiments described above, the indicator .sup.k describes an increasing function of reconstruction noise. In other words, the greater the reconstruction noise, the higher the indicator .sup.k. The optimization algorithm therefore tends to minimize this indicator, in particular on the basis of its gradient .sup.k. The invention may naturally be applied to an indicator describing a decreasing function of reconstruction noise, the indicator decreasing as the reconstruction noise increases. The optimization algorithm then intends to maximize the indicator, in particular on the basis of its gradient. Generally, it is preferable for the noise indicator to follow a monotonic function of the cumulative amplitude of the reconstruction noise in the complex image.

[0142] The invention has been implemented, using the total variation norm, on Chinese hamster ovary (CHO) cells submerged in a CD CHO (Thermo Fisher) culture medium. The sample was placed in a fluidic chamber of 100 m thickness, put at a distance of 8 cm from a light-emitting diode, the spectral band of which was centered on 450 nm. The sample was placed at a distance of 1500 m from a CMOS image sensor of 27483840 pixels. The aperture of the spectral filter 18 has a diameter of 150 m.

[0143] FIG. 5A shows the image I.sub.0 acquired by the image sensor 16. The images of the modulus and of the phase of the reference complex image A.sub.ref, in the plane P.sub.10 of the sample, are shown in FIGS. 5B and 5C, respectively. These images were obtained in 100 iterations. The uniformity of the greyscale levels between each cell attests to the quality of the reconstruction.

[0144] Another example is presented in FIGS. 6A to 6C. In these examples, the sample includes red blood cells diluted in an aqueous solution containing a PBS (phosphate-buffered saline) buffer diluted to 1/400. The sample 10 was placed in a fluidic chamber 15 of 100 m thickness, put at a distance of 8 cm from the light-emitting diode described above, the spectral band of which was centered on 450 nm. The sample was placed at a distance of 1.5 mm from the CMOS image sensor described above. The aperture of the spatial filter 18 was equal to 150 m.

[0145] FIG. 6A shows the image I.sub.0 acquired by the image sensor. The images of the modulus and of the phase of the reconstructed complex image A.sub.z.sup.k=8, in the plane P.sub.10 of the sample, are shown in FIGS. 6B and 6C, respectively. These images were obtained in 8 iterations.

[0146] Whatever the embodiment, it may be necessary to estimate a distance between the detection plane P.sub.0 and the plane P.sub.10 of the sample, in particular when the reference complex image is formed in the plane of the latter. This distance may be known geometrically, or may be estimated by implementing an autofocus algorithm, autofocus algorithms being commonplace in the field of holographic reconstruction.

[0147] The invention will possibly be applied to the observation of a sample by holographic reconstruction, the hologram being obtained either by lens-free imaging, or by defocused imaging. In this case, the hologram is an image acquired by an image sensor, in a plane different from the focal plane of an optical system coupled to the image sensor.

[0148] It will possibly be applied to the characterization of samples in the fields of biotechnology and of diagnostics, but also in the field of food processing, or of the analysis of samples taken from the environment or from industrial processes.