Lensless imaging cell detection method capable of automatically calculating optimal reproduction distance

Abstract

The invention belongs to the technical field of lensless microscopic imaging in the field of bioengineering, and aims to solve the technical problems that in the prior art, only low-concentration cells with dispersed cells can be detected, and the automatic focusing efficiency is relatively low and the accuracy is not high due to manual adjustment of a reproduction distance. A lensless imaging cell detection method capable of automatically calculating an optimal reproduction distance is provided, an obtained original interference image is subjected to numerical reconstruction through an angular spectrum method, length and width information of a to-be-detected sample is determined, a phase unwrapping algorithm is carried out on a digital holographic wrapped phase diagram based on a DCT least square method, and the optimal reproduction distance of the to-be-detected sample is calculated. The obtained phase information is quantitatively converted into thickness information of the sample to be detected, so that the optimal reproduction distance is automatically calculated while reproduction and detection of high-concentration cells are realized, and the accuracy of automatic focusing efficiency is improved.

Description

Technical Field

[0001] The present invention belongs to the field of lensless microscopic imaging technology in the field of bioengineering. Background Art

[0002] Cell detection is of crucial significance in many fields such as biology, medicine, and pharmacy. As an emerging cell detection method, lensless holographic imaging technology has a simple structure, low cost, and does not require staining of samples. Chinese Patent with Publication No. CN106022303A discloses a method for rough classification and counting of freshwater algal cells based on lensless holographic imaging. The holograms of freshwater algae obtained by a lensless holographic imaging device are binarized, and the ratio of the area of all cells to the area of the entire image is calculated. Then, a method based on cell shape features is used to simply classify cell images with dispersed and small numbers of cells.

[0003] However, the above patent can only perform cell detection and simple classification on low-concentration cells with dispersed cells, and most lensless holographic imaging processing methods require manual setting of the reconstruction distance. This is not only cumbersome to operate but also difficult to determine the optimal reconstruction distance, which determines the processing effect and accuracy of the image and affects the quality of subsequent three-dimensional reconstruction and the accuracy of cell counting.

[0004] In summary, there is a need for a lensless holographic imaging cell detection method that can automatically calculate the reconstruction distance to improve the efficiency, accuracy, and reliability of cell detection. For example, Chinese Patent with Publication No. CN118466145A discloses a fast autofocus method for large digital holograms. The reconstructed image amplitude is used as the original matrix, and dimensionality reduction is performed by random projection while retaining features. The reduced-dimensional amplitude is subjected to SVD decomposition, and some singular values are discarded to reduce interference. The L1 norm of the remaining singular values is used as an evaluation function, and the evaluation curve is obtained by traversing the reconstructed image. The clearest position is determined according to the peak value to achieve automatic calculation of the reconstruction distance. However, random projection dimensionality reduction will lose some important information and affect the accuracy of the final focusing. Summary of the Invention

[0005] In order to overcome the technical problems in the above-mentioned prior art that can only perform cell detection on low-concentration cells with dispersed cells, and random projection dimensionality reduction leads to loss of image information, low autofocus efficiency, and low accuracy, the present invention proposes a "lensless imaging cell detection method capable of automatically calculating the optimal reconstruction distance".

[0006] As Figure 1 shown, the lensless imaging cell detection method capable of automatically calculating the optimal reconstruction distance includes the following steps:

[0007] S1: After setting up a lensless digital holographic optical path to perform holographic measurement on the sample to be measured, the interference image generated is received by the image sensor 4, as Figure 4 shown, and is transmitted to the image processing computer as the original interference image, where it is subjected to grayscale processing and converted to double precision. The parameters and the reference light R are set, and the number of horizontal pixels and vertical pixels of the original interference image are recorded.

[0008] R = exp(j·k·(x 0 ·cos(θ 1 ) + y 0 ·cos(θ 2 ))),

[0009] where: exp represents the exponential operation of the natural logarithm e, j is the basic imaginary unit, k = 2π / λ, λ is the central wavelength of the monochromatic light source 3, x 0 is a matrix composed of vectors formed by the lengths of images obtained by dividing the number of vertical pixels of the original interference image by the number of horizontal pixels, and y 0 is a matrix composed of vectors formed by the widths of images obtained by dividing the number of horizontal pixels of the original interference image by the number of vertical pixels. As Figure 3 shown, an x - y - z coordinate system is established with the center of the image sensor 4 as the origin, θ 1 is the angle between the reference light R and the x - axis, and θ 2 is the angle between the reference light R and the y - axis;

[0010] S2: Obtain the optimal reconstruction distance z 0 :

[0011] a. Define the search range and step size of the optimal reconstruction distance z 0 , and determine the number of loops;

[0012] b. Numerically reconstruct the holograms at different reconstruction distances by using the angular spectrum method of Fresnel diffraction integral to obtain a series of reconstructed images;

[0013] c. Unwrap the phase of the obtained reconstructed images in sequence to obtain continuous phase information;

[0014] d. Calculate the standard deviation of the phase gradient and extract the phase of the reconstructed image;

[0015] e. Process the standard deviation of the phase gradient distribution, find the maximum value among the standard deviations of the phase gradients at the recorded different reconstruction distances, and obtain the reconstruction distance z 1 ;

[0016] f. Reset the upper and lower limit values of the reconstruction distance z 1 , obtain a new step size, and repeat steps b - f to obtain the optimal reconstruction distance z 0 ;

[0017] S3: Determine the optimal reproduction distance z 0 After that, the angular spectrum method is used to reproduce the original interference image obtained in S1. As Figure 5 shown, the length and width of the sample to be measured are determined through reproduction, and the complex amplitude distribution of the object light field representing the light intensity information of the sample to be measured is obtained:

[0018]

[0019] wherein, is the complex amplitude distribution of the object light field, λ is the central wavelength of the monochromatic light source 3, F -1 is the inverse Fourier transform, F is the Fourier transform, (x, y) is the diffraction plane coordinate, U 0 (x, y) is the complex amplitude of the object plane, exp is the exponential operation of the natural logarithm e, i is the imaginary unit, (f xi , f yi ) is the frequency domain coordinate, z 0 is the optimal reproduction distance;

[0020] S4: Calculate the wrapped phase Ψ(x, y) of the corresponding sample to be measured through the phase extraction function. The specific formula is as follows:

[0021]

[0022] In the formula, is the complex amplitude distribution of the object light field. The phase information in the object light field is obtained through the arctangent function, where represents the real part of the object light field,

[0023]

[0024] Ψ 0 (x, y) = Ψ(x, y) + 2k(x, y)π,

[0024] In the formula: Ψ(x, y), Ψ 0 (x, y) respectively represent the wrapped phase and the original phase of the sample to be measured. The value of k(x, y) is an integer,

[0025]

[0026] In the formula: Ψ 0 (x, y) is the original phase of the sample to be measured, Δl is the optical path difference after the monochromatic light source passes through the sample to be measured and the solution medium, λ is the central wavelength of the monochromatic light source 3,

[0027] Δl = d(x, y) × (n p-n medium )

[0028] Subsequently, the thickness of the sample to be measured is obtained Wherein:

[0029] d(x,y) represents the thickness of the sample to be measured, and n p represents the refractive index of the sample to be measured, and n medium is the refractive index of the solution medium where the sample to be measured is located, and Ψ 0 (x,y) represents the original phase of the object to be measured, and λ is the central wavelength of the monochromatic light source 3

[0030] Technical effects:

[0031] In order to overcome the problem that only low-concentration cells with cell dispersion can be detected in the prior art, the present invention determines the length and width information of the sample to be measured through angular spectrum method cell reproduction numerical reconstruction, and performs a phase unwrapping algorithm on the digital holographic wrapped phase diagram based on the least squares method of DCT. As Figure 10 shown, it is a schematic diagram of the three-dimensional phase unwrapping effect of onion epidermal cells, without calculating the residual points or quality map of the wrapped phase diagram, and quantitatively converting the obtained phase information into the thickness information of the sample to be measured by using the optical path difference and phase conversion in optics; S2 in the present invention solves the technical problems of image information loss caused by random projection dimensionality reduction, low efficiency and low accuracy of manual adjustment of the reproduction distance for autofocusing, and determines the optimal reproduction distance z by using the gradient standard deviation of the calculated phase distribution 0 , and the judgment method that the greater the gradient standard deviation, the better the reproduction effect is simple and effective. Compared with the prior art, this method not only has a simple and practical algorithm, but also can more accurately highlight the sample to be measured, quickly determine the optimal reproduction position, and improve the accuracy of the autofocusing efficiency

[0032] As Figure 5 shown, it is a reproduction diagram reconstructed by the angular spectrum method after the sample to be measured finds the optimal distance through the automatic calculation of the reproduction distance algorithm. Compared with Figure 7 , the accuracy of the optimal reproduction distance algorithm proposed by the present invention can be verified. As Figure 8 and Figure 9 shown, they are the interference image of onion epidermal cells and the reproduction diagram reconstructed by the angular spectrum method after the optimal reproduction distance algorithm, which not only verifies that the present invention can realize the reproduction and detection of high-concentration cells, but also verifies the accuracy of the proposed optimal reproduction distance algorithm Brief description of the drawings

[0033] Figure 1 This is the overall flow block diagram of the present invention and is used as the abstract drawing

[0034] Figure 2 It is a schematic diagram when the original interference image is obtained by the lensless digital holographic optical path

[0035] Figure 3 The angle θ formed by the reference light R 1 and θ 2 Schematic diagram of the geometric meaning.

[0036] Figure 4 Schematic diagram of the original interference image of the sample to be measured obtained under the lensless digital holographic optical path.

[0037] Figure 5 The angular spectrum method reproduction diagram of the original interference image of the present invention.

[0038] Figure 6 The angular spectrum method reproduction diagram of the filtered image of the present invention.

[0039] Figure 7 The angular spectrum method reproduction diagram of the original interference image with the reproduction distance manually set.

[0040] Figure 8 Schematic diagram of the interference image effect of onion epidermal cells in the embodiment.

[0041] Figure 9 Schematic diagram of the interference image reproduction effect of onion epidermal cells in the embodiment.

[0042] Figure 10 Schematic diagram of the three-dimensional phase unwrapping effect of onion epidermal cells in the embodiment.

[0043] Figure 11 Schematic diagram of the phase-thickness conversion principle. Detailed implementation manner

[0044] For the lensless digital holographic optical path, as Figure 2 shown, the sample to be measured is placed between two cover glasses 1 to form a thin layer 2, and the light emitted by the monochromatic light source 3 is vertically incident on the thin layer 2, and the generated interference image is received by the image sensor 4 and transmitted to the image processing computer as the original interference image.

[0045] The monochromatic light source 3 is a monochromatic LED with a central wavelength of 630.0 nm. Its power supply voltage is controlled to be 3 V and the current is 0.1 A to make the light source vertically irradiate the thin layer 2. A micropore can be set below the monochromatic light source 3 to limit the light wave, and the diameter of the micropore is reduced to increase the spatial coherence of the light source. The diameter of the micropore is set to 100 μm. The image sensor 4 is a CMOS sensor with a pixel size of 1.4 μm × 1.4 μm; the specific parameters are the central wavelength λ of the monochromatic light source 3 and the distance z from the thin layer 2 to the image sensor 4 2 .

[0046] Such as Figure 3As shown in the figure, an x-y-z coordinate system is established with the center of the image sensor 4 as the origin, and θ 1 is the angle between the reference light R and the x-axis, and θ 2 is the angle between the reference light R and the y-axis; when the Z-axis represents the reference light R, θ 1 and θ 2 are preferably 90°.

[0047] Since the optimal reproduction distance is extremely close to the distance from the sample to be measured to the image sensor, and there is a protective glass about 0.5 mm thick on the image sensor, the optimal reproduction distance is between 0.5 mm and 1 mm. The specific expression of the number of cycles times is:

[0048]

[0049] In the formula: z 0max is the upper limit value 1 mm of the reproduction distance, z 0min is the lower limit value 0.5 mm of the reproduction distance, is rounding up, z 0n is the step size. First, calculate how many z 0n are included in the search range, then round up and add 1 to get the number of cycles. Inside the cycle, numerically reconstruct the interference images at different reproduction distances through the angular spectrum method to obtain a series of reproduced images. Unwrap the phase of the obtained reproduced images in turn to obtain continuous phase information, and increase the accuracy of the algorithm for automatically calculating the optimal reproduction distance.

[0050] Calculate the standard deviation of the phase gradient. The gradient refers to the directional derivative of a certain function at that point along the direction where the maximum value is obtained. This direction is the gradient direction, and this value is the modulus of the gradient. The specific calculation method is:

[0051]

[0052] is the modulus of the gradient, G x (r, c) is the gradient component along the x direction at the coordinate (r, c), G y (r, c) is the gradient component along the y direction at the coordinate (r, c). Then extract the phase of the reproduced image, and then process the standard deviation of the gradient of the phase distribution through the algorithm. Find the maximum value among the standard deviations of the gradient amplitudes at different recorded reproduction distances to obtain the reproduction distance z 1 , then reset the upper limit value and the lower limit value of the reproduction distance, obtain a new step size, and repeat the above operations to determine the optimal reproduction distance z 0 .

[0053] Before the original interference image obtained in S1 is reproduced using the angular spectrum method in S3, filtering processing can be performed. Further, an ideal high-pass filter is used to filter the processed original interference image to remove the influence of background pixel values and perform spectral component adjustment, as Figure 6 shown.

[0054] The ideal high-pass filter G 1 (u, v) has the following specific expression:

[0055]

[0056] In the formula, u and v are the horizontal and vertical pixel numbers of the original interference image, is the distance from the point (u, v) to the origin, D 0 is the distance from the cut-off frequency point to the origin. In the ideal high-pass filter, the value of D 0 is 0.0015 times the maximum frequency, and the value of the maximum frequency is equal to the value of the image diagonal.

[0057] In S3, the angular spectrum method is used to reproduce the original interference image obtained in S1, as Figure 5 shown. Based on the principle of light diffraction, numerical reconstruction is carried out. The interference fringe information recorded by the image sensor is transmitted to the image processing computer, and the complex amplitude distribution of the light wave in the reconstructed image plane is obtained through numerical simulation to restore the true morphology of the sample to be measured. As an accurate expression of the scalar diffraction theory, the angular spectrum method of Fresnel diffraction can achieve accurate calculation of the diffraction light field in free space. The angular spectrum method is expressed as:

[0058]

[0059] Among them, is the complex amplitude distribution of the object light field, λ is the central wavelength of the monochromatic light source 3, F -1 is the inverse Fourier transform, F is the Fourier transform, (x, y) is the diffraction plane coordinate, U 0 (x, y) is the complex amplitude of the object plane, exp is the exponential operation of the natural logarithm e, i is the imaginary unit, (f xi , f yi ) is the frequency domain coordinate, z 0 is the optimal reconstruction distance.

[0060] In the numerical reconstruction of holograms, the angular spectrum method has a total of two Fourier transforms, with a relatively large amount of calculation and computational cost. However, the algorithms for automatically calculating the optimal reconstruction distance in the present invention combined with Fresnel and convolution methods are only applicable to holograms of suitable sizes. Only the reconstruction combined with the angular spectrum method is not affected by the reconstruction size, and the obtained optimal reconstruction distance can focus on the reconstruction of cell holograms under different reconstruction sizes.

[0061] S4. Calculate the wrapped phase Ψ(x, y) of the corresponding sample to be measured through the phase extraction function. The specific formula is as follows: wherein,

[0062]

[0063] In the formula, is the complex amplitude distribution of the object light field. The phase information in the object light field is obtained through the arctangent function, where represents the real part of the object light field, represents the imaginary part of the object light field. Then, the phase unwrapping algorithm is used to obtain the original phase of the sample to be measured. According to the periodic property of the arctangent function, the variation range of the phase information Ψ(x, y) is between (-π / 2, π / 2). Since U 0 (x, y) is in complex form and its variation range exists in four quadrants, the variation range of Ψ(x, y) is extended from (-π / 2, π / 2) to (-π, π). For a sample to be measured with a certain thickness, if its original phase value is outside the range of (-π, π), the phase value directly obtained from the reconstructed object light field will still fall within the range of (-π, π). Therefore, for any phase value Ψ(x, y), a ±2kπ compensation is required using the phase unwrapping algorithm:

[0064] Ψ 0 (x, y) = Ψ(x, y) + 2k(x, y)π,

[0065] In the formula: Ψ(x, y) and Ψ 0 (x, y) represent the wrapped phase and the original phase of the sample to be measured respectively. The value of k(x, y) is an integer. If the original phase of the sample to be measured is to be obtained, the wrapped phase obtained needs to be unwrapped. This process is called phase unwrapping, and the algorithm used is the phase unwrapping algorithm.

[0066] During the recording process of digital holography, the optical path differences between the object light wave and the reference light wave caused by different refractive indexes inside the sample to be measured are also different. The obtained phase value is closely related to the true thickness of the object and the refractive index of its internal structure. As Figure 11 shown, using the knowledge of the conversion between optical path difference and phase in optics, the obtained phase information is quantitatively converted into the thickness information of the sample. The phase reflects the relationship between the refractive indexes of the internal structure of the sample to be measured. Therefore, the topography formed on the surface of the obtained phase information is consistent with the shape formed by integrating the refractive index of the sample to be measured along the direction of the reference light wave inside the sample.

[0067] The relationship between the optical path difference and the phase value in optical principles can be expressed as:

[0068]

[0069] In the formula: Ψ 0(x,y) is the original phase of the sample to be measured, Δl is the optical path difference after the monochromatic light source passes through the sample to be measured and the solution medium, λ is the central wavelength of the monochromatic light source 3. According to the conversion relationship between refractive index, light speed and wavelength, the optical path difference can be expressed as

[0070] Δl = d(x,y)×(n p -n medium ),

[0071] Then the thickness of the sample to be measured is obtained:

[0072] After automatically calculating the optimal reproduction distance of the present invention, taking the cell detection of onion epidermis cells as an example, combined with the attached Figures 8 - 10 The corresponding effect schematic diagram is given. The above is only an example clearly illustrating the present invention, rather than a limitation on the implementation mode. For those of ordinary skill in the art, other different forms of changes can be made based on the above description. It is not necessary and impossible to enumerate all implementation modes here, and the obvious changes or variations derived therefrom are still within the protection scope of the present invention.

Claims

1. A lens-free imaging cell detection method capable of automatically calculating an optimal reproduction distance, comprising the following steps: S1: The interference image generated after the lensless digital holographic optical path is constructed and the holographic measurement of the sample to be measured is received by the image sensor (4), and transmitted to the image processing computer as the original interference image, grayscale processing is performed and converted into double precision, parameters and reference light R are set, and the horizontal pixel number and vertical pixel number of the original interference image are recorded. R=exp(j·k·(x0·cos(θ1)+y0·cos(θ2))), In the formula: exp represents the exponential operation of the natural logarithm e, j is the basic imaginary unit, k=2π / λ, λ is the central wavelength of the monochromatic light source (3), x0 is a matrix composed of vectors composed of the length of the image divided equally by the number of horizontal pixels of the original interference image, y0 is a matrix composed of vectors composed of the width of the image divided equally by the number of horizontal pixels of the original interference image, an xyz coordinate system is established with the center of the image sensor (4) as the origin, θ1 is the angle between the reference light R and the x-axis, and θ2 is the angle between the reference light R and the y-axis; it is characterized in that, S2: Get the best reproduction distance z0: a. Define the search range and step size of the optimal reproduction distance z0 and determine the number of cycles; b. numerically reconstructing the holograms at different reconstruction distances by using the angular spectrum method of Fresnel diffraction integral to obtain a series of reconstruction images; c. performing phase unwrapping on the obtained reconstructed images in sequence to obtain continuous phase information; d. Calculate the phase gradient standard deviation and extract the phase of the reconstructed image; e. Process the gradient standard deviation of the phase distribution, find the maximum value among the phase gradient standard deviations recorded at different reproduction distances, and obtain the reproduction distance z1; f. Reset the upper and lower limits of the reproduction distance z1 to obtain a new step length, and repeat steps b to f to obtain the optimal reproduction distance z0; S3: After determining the optimal reproduction distance z0, the original interference image obtained in S1 is reproduced using the angular spectrum method. The length and width of the sample to be tested are determined by reproduction, and the complex amplitude distribution of the object light field representing the light intensity information of the sample to be tested is obtained: in, is the complex amplitude distribution of the object light field, λ is the central wavelength of the monochromatic light source (3), F -1 is the inverse Fourier transform, F is the Fourier transform, (x, y) is the diffraction surface coordinate, U0(x, y) is the complex amplitude of the object surface, exp is the exponential operation of the natural logarithm e, i is the imaginary unit, (f xi ,f yi ) is the frequency domain coordinate, z0 is the reproduction distance; S4: Calculated by phase extraction function The corresponding wrapped phase Ψ(x,y) of the sample to be tested is as follows: In the formula, is the complex amplitude distribution of the object light field. The phase information in the object light field is obtained by the inverse tangent function, where represents the real part of the object light field, Represents the imaginary part of the object light field, and then uses the phase unwrapping algorithm to obtain the original phase of the sample to be tested. For any wrapped phase Ψ(x, y), the unwrapping algorithm needs to compensate ±2kπ: Ψ0(x,y)=Ψ(x,y)+2k(x,y)π, Where: Ψ(x, y) and Ψ0(x, y) represent the wrapped phase and original phase of the sample to be tested respectively, and the value of k(x, y) is an integer. Where: Ψ0(x, y) is the original phase of the sample to be tested, Δl is the optical path difference of the monochromatic light source after passing through the sample to be tested and the solution medium, λ is the central wavelength of the monochromatic light source (3), Δl=d(x,y)×(n p -n medium ), Then the thickness of the sample to be tested is obtained in: d(x,y) represents the thickness of the sample to be tested, n p Represents the refractive index of the sample to be tested, n medium is the refractive index of the solution medium in which the sample to be tested is located, Ψ0(x, y) represents the original phase of the object to be tested, and λ is the central wavelength of the monochromatic light source (3).

2. The lensless imaging cell detection method capable of automatically calculating the optimal reproduction distance according to claim 1, characterized in that: The lensless digital holographic optical path is specifically as follows: a sample to be tested is placed between two cover glasses (1) to form a thin layer (2); light emitted by a monochromatic light source (3) vertically illuminates the thin layer (2); a generated interference image is received by an image sensor (4) and transmitted to an image processing computer as an original interference image; the parameters are specifically set as the central wavelength λ of the monochromatic light source (3) and the distance z2 from the thin layer (2) to the image sensor (4).

3. The lensless imaging cell detection method capable of automatically calculating the optimal reproduction distance according to claim 2, characterized in that: The monochromatic light source (3) is a monochromatic LED with a central wavelength of 630.0 nm. The power supply voltage is controlled to be 3 V and the current is controlled to be 0.1 A so that the light emitted by the monochromatic light source (3) vertically irradiates the thin layer (2). A microhole is arranged below the monochromatic light source (3). The diameter of the microhole is set to be 100 μm. The image sensor (4) is a CMOS sensor with a pixel size of 1.4 μm×1.4 μm.

4. The lensless imaging cell detection method capable of automatically calculating the optimal reproduction distance according to claim 1, characterized in that: The specific expression of the number of cycles in S2 is: Where: times is the number of cycles, z 0max The upper limit of the reproduction distance is 1mm, z 0min The lower limit of the reproduction distance is 0.5 mm. To round up, z 0n is the step length.

5. The lensless imaging cell detection method capable of automatically calculating the optimal reproduction distance according to claim 1, characterized in that: The specific expression for calculating the phase gradient standard deviation in S2 is: Where: is the modulus of the gradient, G x (r,c) is the gradient component along the x direction at the coordinate (r,c), G y (r,c) is the gradient component along the y direction at coordinate (r,c).

6. The lensless imaging cell detection method capable of automatically calculating the optimal reproduction distance according to claim 1, characterized in that: It also includes filtering the original interference image obtained by S1 before reproducing it using the angular spectrum method. The specific expression of the ideal high-pass filter G1(u,v) is: Where u and v are the horizontal and vertical pixel numbers of the original interference image. is the distance from the point (u, v) to the origin, D0 is the distance from the cutoff frequency point to the origin, and the value of D0 in an ideal high-pass filter is 15 thousandths of the maximum frequency, and the value of the maximum frequency is equal to the value of the image diagonal.

Images