An adaptive aberration compensation quantitative phase imaging method based on annular illumination

Abstract

The application discloses a kind of adaptive aberration compensation quantitative phase imaging methods based on annular illumination, through to the sample to be measured 6 to 12 light intensity maps of illumination angle equal to objective numerical aperture are shot, iteration reconstruction is carried out, and the phase of sample to be measured and system aberration are solved simultaneously finally, real-time adaptive aberration compensation quantitative phase imaging of sample to be measured is realized.The application can solve the problem of focal length drift caused by environmental temperature change, mechanical structure instability and other reasons when observing living cells for a long time, and real-time correction of time-varying aberration occurring in the imaging process.

Description

Technical Field

[0001] The present invention belongs to optical microscopy imaging technology, and in particular to a self-adaptive aberration compensation quantitative phase imaging method based on annular illumination. Background Art

[0002] Quantitative phase imaging (QPI) is a technique that measures the phase of an object based on the intensity of the imaged light. It holds great promise in the biomedical field. Over the past decade, numerous QPI techniques have made significant progress, including digital holography, the intensity transfer equation, differential phase contrast, and Fourier stack imaging. Unlike traditional fluorescence imaging, QPI is label-free and stain-free, making it ideal for real-time observation of living cells. However, it is important to note that when using QPI to measure live cells over extended periods, environmental disturbances such as room temperature fluctuations, undulations in the culture medium surface, and vertical mechanical displacement of the stage can significantly impact measurement results. These factors can cause real-time changes in the imaging system's aberrations, ultimately compromising the quality and accuracy of phase imaging. Therefore, adaptive aberration correction is essential for maintaining high-quality, long-term phase imaging of living cells.

[0003] Adaptive optics (AO) is a common method used in microscopy to remove the effects of aberrations and achieve high-resolution imaging. Aberration compensation based on adaptive optics is primarily divided into two steps. First, the aberrations are actively measured, such as using a Shakman wavefront sensor, interferometry, phase retrieval, and point spread function measurement. The second step is to compensate for the measured aberrations. A common approach in adaptive optics is to use a deformable mirror (DM) and a spatial light modulator (SLM). However, all of the above adaptive optics-based methods require the introduction of additional hardware into the system, increasing system complexity and cost. In contrast, Fourier stacking (FPM), a computational imaging technique, can exploit the redundancy of the acquired data to computationally solve the system aberrations without the need for any additional hardware. However, traditional FPM methods require a large amount of additional data to simultaneously resolve information such as the object's phase and aberrations, resulting in slow imaging speeds that cannot meet the needs of real-time observation of living cells. Summary of the Invention

[0004] The purpose of the present invention is to propose an adaptive aberration compensation quantitative phase imaging method based on annular illumination, which solves the problem of focal drift caused by environmental temperature changes, mechanical structure instability, etc. during long-term observation of living cells, and can correct the time-varying aberrations that occur during the imaging process in real time.

[0005] The technical solution for achieving the purpose of the present invention is: a quantitative phase imaging method based on adaptive aberration compensation and annular illumination, which specifically comprises the following steps:

[0006] Step 1: Use a ring LED as the lighting source and calibrate the LED brightness;

[0007] Step 2: Light up 6 to 12 LEDs distributed in a circular pattern with equal spacing. Collect a series of original light intensity images of the sample under test at different lighting angles under the synchronous trigger signal of the camera.

[0008] Step 3: Preprocess the original light intensity image, including threshold denoising and image brightness correction;

[0009] Step 4: Initialize the spectrum of the sample to be tested;

[0010] Step 5: Use the normalized light intensity map to iteratively update the pupil function and sample spectrum in turn. When the difference between the cost functions of two iterations is less than the set threshold, stop the iteration and use the current sample spectrum as the final calculation result.

[0011] Preferably, the specific method for calibrating the LED brightness is:

[0012] Adjust the LED to a fixed height, light up 6 to 12 LEDs with equal spacing in a ring, and collect a series of images with different lighting angles under the camera synchronization trigger signal. B 1(x,y),I B 2(x,y),...,I B N (x, y), calculate the brightness normalization coefficient of each LED, specifically: for images I at different lighting angles B 1(x,y),I B 2(x,y),...,I B N (x, y) is averaged and normalized to obtain the brightness normalization coefficients α1, α2, ..., α N .

[0013] Preferably, the specific method for preprocessing the original light intensity image is:

[0014] According to the average value of the camera dark current noise as the threshold, the original image is threshold denoised, and the original light intensity map of the sample I1(x,y),I2(x,y),...,I N (x,y) are divided by the brightness normalization coefficient.

[0015] Preferably, the specific method for initializing the spectrum of the sample to be measured is:

[0016] Normalized light intensity images of samples I1(x,y), I2(x,y),...,I N (x,y) is Fourier transformed to obtain the corresponding spectrum S1(u,v), S2(u,v),..., S N (u,v), filter the spectra separately and add them together to get the initialized sample spectrum U 0 (u,v), specifically:

[0017]

[0018] Where N is the number of light intensity images, f n (u,v) is the filter.

[0019] Preferably, the filter for filtering the spectrum is:

[0020]

[0021] Where θ n is the angle of the nth LED relative to the x-axis of the camera pixel plane, (u, v) is the spectrum coordinate

[0022] coordinate system.

[0023] Preferably, the pupil function and the sample spectrum are iteratively updated in sequence using the normalized light intensity map.

[0024] The steps are:

[0025] Update the spectrum U(u,v), the specific formula is:

[0026] U j+1 (uu n ,vv n )=U j (uu n ,vv n )+conj(P j (u,v)) / |P j (u,v) 2 ·[U n j,update (u,v)-U j n (u,v)]

[0027] Where j represents the number of iterations, n represents the sequence number of the updated light intensity map and the corresponding LED lighting angle under the current iteration, (u n ,v n ) represents the coordinates of the spectrum displacement corresponding to the nth light intensity image, U n (u, v) represents the spectrum within the aperture corresponding to the nth illumination angle, specifically:

[0028] U j n (u,v)=U j n (uu n ,vv n )P j (u,v)

[0029] For the light intensity diagram I n The spectrum of (x,y) after spatial constraints:

[0030]

[0031] Among them, u j n (x,y) and U j n (u,v) are Fourier transform pairs;

[0032] Pupil function P j (u,v) is updated by the system aberration A in the iterative update j (u,v) and the low-pass filter circ(NA / λ) are expressed as:

[0033] P j (u,v)=circ(NA / λ)exp(iA j (u,v))

[0034] Where, the system aberration A j (u,v) is calculated from the first 15 Zernike polynomials i The weighted sum of (u,v) represents p j i is the weighted coefficient of the polynomial, and i represents the serial number of the Zernike polynomial coefficient:

[0035]

[0036] By p j i Update the pupil function P(u,v) and initialize the polynomial coefficients p when j=0. 0 i = 0, in the jth iteration, the polynomial coefficient p j i The update formula is:

[0037]

[0038] Preferably, the cost function is specifically:

[0039]

[0040] Compared with the prior art, the present application has the following advantages: (1) Compared with the traditional quantitative phase imaging technology, the present application adds the function of real-time aberration compensation without adding additional hardware in the imaging system. (2) The present application only needs to shoot 6-12 light intensity images to realize quantitative phase imaging and aberration correction, has high imaging efficiency, and can be applied to real-time dynamic imaging of living cells.

[0041] The present application will be described in further detail below with reference to the accompanying drawings. BRIEF DESCRIPTION OF DRAWINGS

[0042] Figure 1 is an illumination mode diagram using different numbers of annular LEDs. The white dots in the diagram represent the spatial frequency positions corresponding to the LEDs in the frequency domain, and the white dashed line represents the numerical aperture NA of the objective lens.

[0043] Figure 2 is a flowchart of the adaptive aberration compensation quantitative phase imaging method based on annular illumination.

[0044] Figure 3 is a diagram of the aberration-corrected quantitative phase imaging result and the recovered imaging system aberration achieved by the method. The small diagrams (a2) (a3) and (b2) (b3) on the right represent the cell phase images obtained after aberration correction and without aberration correction.

[0045] Figure 4 is a long-time imaging result of HeLa living cells for 48 hours using the method and the recovered time-varying aberration. DETAILED DESCRIPTION

[0046] The present application will be described in further detail below with reference to the accompanying drawings.

[0047] The present application is an adaptive aberration compensation quantitative phase imaging method based on annular illumination, comprising the following steps:

[0048] Step 1: LED brightness calibration: using annular LEDs as the illumination light source, adjusting the LEDs to a specific height and fixing them, sequentially lighting 6 to 12 annular LEDs distributed at equal intervals, collecting an image for each lighting, and collecting a series of images I B 1(x,y), I B 2(x,y),..., I B N (x,y) to calculate the brightness normalization coefficients α1, α2,..., α N In the imaging process, annular LEDs with a radius of r are used as the illumination light source, and the numerical aperture of the objective lens is NA objWhen imaging with the objective lens, the LED needs to be adjusted to a specific height h to meet Figure 1 Figure 2 shows the lighting patterns using different numbers of ring-shaped LEDs.

[0049] Step 2: Original image acquisition: Use a ring LED as the illumination source, light up 6 to 12 ring-shaped LEDs with equal spacing in sequence, and collect a series of original light intensity images I1(x,y), I2(x,y), ..., I at different illumination angles of the sample under test under the synchronous trigger signal of the camera. N (x,y).

[0050] Step 3: Preprocessing the original image: The original image preprocessing includes threshold denoising and image brightness correction. First, the average value of the camera dark current noise is used as the threshold. The pixel values ​​of the original image are compared with the threshold. The pixel values ​​less than the threshold are set to zero to complete the threshold denoising. Then, according to the brightness normalization coefficients α1, α2, ..., α obtained in step 1, the original image is denoised. N , the original light intensity map of the sample to be tested I1(x,y),I2(x,y),...,I N (x, y) is normalized for brightness. The specific operation is as follows: calculate the image I of different lighting angles under the blank field of view respectively. B 1(x,y),I B 2(x,y),...,I B N The average value of (x, y) is then normalized to obtain the brightness normalization coefficients α1, α2, ..., α N Then the original light intensity map of the sample I1(x,y),I2(x,y),...,I N (x,y) are divided by the brightness normalization coefficient respectively.

[0051] Step 4: Initialize the spectrum of the sample to be tested: normalize the light intensity graphs I1(x,y), I2(x,y0,...,I N (x,y) is Fourier transformed to obtain the corresponding spectrum S1(u,v), S2(u,v),..., S N (u,v), filter the spectra separately and add them together to get the initialized sample spectrum U 0 (u,v). Filter f for each original light intensity map spectrum n The calculation formula for (u,v) is as follows:

[0052]

[0053] Among them, θ n is the angle between the nth LED and the x-axis of the camera pixel plane. The calculation formula for the initialization sample spectrum U0 is as follows:

[0054]

[0055] Step 5: Iterative reconstruction: In each round of iterative reconstruction, the normalized intensity map I1(x,y),I2(x,y),...,I N (x,y) updates the pupil function P(u,v) and sample spectrum U(u,v) in turn.

[0056] The update formula for the spectrum U(u,v) in iterative reconstruction is:

[0057] U j+1 (uu n ,vv n )=U j (uu n ,vv n )+conj(P j (u,v)) / |P j (u,v) 2 ·[U n j,update (u,v)-U j n (u,v)]

[0058] Among them, j represents the iteration round number, n represents the sequence number of the updated light intensity map and the corresponding LED lighting angle under the current iteration round number, (u n ,v n ) represents the coordinates of the spectrum displacement corresponding to the nth light intensity image. n (u,v) represents the spectrum within the aperture corresponding to the nth illumination angle:

[0059] U j n (u,v)=U j n (uu l ,vv l )P j (u,v)

[0060] For the light intensity diagram I n The spectrum of (x,y) after spatial constraints:

[0061]

[0062] Among them, u j n (x,y) and U j n (u,v) are Fourier transform pairs.

[0063] Pupil function P j(u,v) is updated by the system aberration A in the iterative update j (u,v) and the low-pass filter circ(NA / λ) are expressed as:

[0064] P j (u,v)=circ(NA / λ)exp(iA j (u,v))

[0065] Where, the system aberration A j (u,v) is calculated from the first 15 Zernike polynomials i (u,v) [6] The weighted sum of j i is the weighted coefficient of the polynomial, and i represents the serial number of the Zernike polynomial coefficient:

[0066]

[0067] By p j i Update the pupil function P(u,v) when j=0, initialize the polynomial coefficient p 0 i = 0. In the jth iteration, the polynomial coefficient p j i The update formula is:

[0068]

[0069] When the pupil function P(u,v) and the sample spectrum U(u,v) are updated, the jth iteration ends, and then the cost function cost is calculated as the convergence criterion. j+1 :

[0070]

[0071] When the cost function cost j Satisfaction (cost j -cost j+1 ) / cost j When <0.01, the iteration stops at the jth round, and the sample spectrum U j For the final calculation result, the sample phase is Φ(x,y)=angle(F -1 (U j (u,v))); otherwise, continue with the j+1th iteration.

[0072] The present invention is based on Fourier stack imaging technology and uses only the light intensity images captured by a small number of LEDs distributed in a ring with an illumination angle equal to the numerical aperture of the objective lens as raw data to achieve quantitative phase imaging and aberration correction. It can realize dynamic aberration compensation during real-time imaging of living cells, and meet the requirements of long-term stable high-quality imaging of living cells.

[0073] Figure 2 is a schematic diagram of the workflow of this method.

[0074] Figure 3 Comparison of the results before and after using this method to correct the aberration of the imaging system. Figures (a1) and (b1) show the results of quantitative phase imaging of living cells without aberration correction and the reconstructed aberration. Figures (a2) and (b2) are magnified images of the local area, and Figures (a3) ​​and (b3) are comparisons of the phase imaging results without aberration correction. The figures show that the proposed method can effectively improve the imaging quality of living cells.

[0075] Figure 4 This figure shows the results of real-time aberration compensation for long-term live cell imaging using this method. The figure shows the cell phase and system aberration at different times. It can be seen that thanks to the real-time aberration correction, the cell phase maintains a high imaging quality throughout the long-term imaging process.

Claims

1. A quantitative phase imaging method with adaptive aberration compensation based on annular illumination, characterized in that: The specific steps are: Step 1: Use a ring LED as the lighting source and calibrate the LED brightness; Step 2: Light up 6 to 12 LEDs distributed in a circular pattern with equal spacing. Collect a series of original light intensity images of the sample under test at different lighting angles under the synchronous trigger signal of the camera. Step 3: Preprocess the original light intensity image, including threshold denoising and image brightness correction; Step 4: Initialize the spectrum of the sample to be tested; Step 5: Use the normalized light intensity map to iteratively update the pupil function and sample spectrum. When the difference between the cost functions of two iterations is less than the set threshold, stop the iteration and use the current sample spectrum as the final calculation result. The specific steps are as follows: Update the spectrum U(u,v), the specific formula is: U j+1 (u-u n ,v-v n )=U j (u-u n ,v-v n )+conj(P j (u,v)) / |P j (u,v)| 2 ·[U n j,update (u,v)-U j n (u,v)] Where j represents the number of iterations, n represents the serial number of the updated light intensity map and the corresponding LED lighting angle under the current iteration, (u n ,v n ) represents the coordinates of the spectrum displacement corresponding to the nth light intensity image, U n (u, v) represents the spectrum within the aperture corresponding to the nth illumination angle, specifically: U j n (u,v)=U j n (uu n ,vv n )P j (u,v) U j n update (u,v) is the light intensity image I n The spectrum of (x,y) after spatial constraints: Among them, u j n (x,y) and U j n (u,v) are Fourier transform pairs; Pupil function P j (u,v) is updated by the system aberration A in the iterative update j (u,v) and low-pass filter circ(NA / λ) means: P.S j (u,v)circ(NA / λ)exp(iA j (u,v)) Where, the system aberration A j (u,v) is calculated from the first 15 Zernike polynomials i The weighted sum of (u,v) represents p j i is the weighted coefficient of the polynomial, and i represents the serial number of the Zernike polynomial coefficient: By p j i Update the pupil function P(u,v) and initialize the polynomial coefficients p when j=0. 0 i = 0, in the jth iteration, the polynomial coefficient p j i The update formula is: The cost function is specifically:

2. The method for quantitative phase imaging with adaptive aberration compensation based on annular illumination according to claim 1, wherein: The specific method for calibrating LED brightness is: Adjust the LED to a fixed height, light up 6 to 12 LEDs with equal spacing in a ring, and collect a series of images with different lighting angles under the camera synchronization trigger signal. B 1(x,y),I B 2(x,y),...,I B N (x, y), calculate the brightness normalization coefficient of each LED, specifically: for images I at different lighting angles B 1(x,y),I B 2(x,y),...,I B N (x, y) is averaged and normalized to obtain the brightness normalization coefficients α1, α2, ..., α N .

3. The method for quantitative phase imaging with adaptive aberration compensation based on annular illumination according to claim 1, wherein: The specific method for preprocessing the original light intensity map is: According to the average value of the camera dark current noise as the threshold, the original image is threshold denoised, and the original light intensity map of the sample I1(x,y),I2(x,y),...,I N (x,y) are divided by the brightness normalization coefficient.

4. The method for quantitative phase imaging with adaptive aberration compensation based on annular illumination according to claim 1, wherein: The specific method for initializing the spectrum of the sample to be tested is: Normalized light intensity images of samples I1(x,y), I2(x,y),...,I N (x,y) is Fourier transformed to obtain the corresponding spectrum S1(u,v), S2(u,v),..., S N (u,v), filter the spectra separately and add them together to get the initialized sample spectrum U 0 (u,v), specifically: Where N is the number of light intensity images, f n (u,v) is the filter.

5. The method for quantitative phase imaging with adaptive aberration compensation based on annular illumination according to claim 4, characterized in that: The filter used to filter the spectrum is: Where θ n is the angle of the nth LED relative to the x-axis of the camera pixel plane, and (u, v) is the coordinate in the spectrum coordinate system.

Images