Fast 3D deconvolution microscopy imaging method based on mismatched backprojection

Abstract

The present invention proposes a fast three-dimensional deconvolution microscopy imaging method based on mismatched back projection. The method mainly includes: the method uses deconvolution to restore the image in a fluorescence incoherent mode, and uses a Wiener-Butterworth filter to generate an unmatched back projector. Different from the back projector that directly uses a mirror image to generate a point spread function (PSF) in the traditional RL algorithm, the method can achieve three-dimensional image reconstruction with high precision and high efficiency, while restoring the image to the theoretical resolution limit, greatly improving the imaging speed, and reducing artifacts caused by the increase in the number of deconvolutions; CUDA is used to parallelize the three-dimensional deconvolution imaging method, and the speed is increased by 30-40 times again through the powerful floating-point computing power of the GPU, so as to achieve efficient and fast calculation of large-scale data. The fast three-dimensional deconvolution microscopy imaging method provided by the present invention supports real-time three-dimensional deconvolution and can meet the urgent needs of fields such as dynamic life science research and medical imaging, and has broad application prospects.

Description

Technical Field

[0001] The invention belongs to the technical field of optical microscopic imaging, and in particular relates to a fast three-dimensional deconvolution microscopic imaging method based on mismatched back projection. Background Art

[0002] With the increasingly urgent needs of life sciences, optical microscopy has become an indispensable and important tool to promote life science research and modern clinical disease diagnosis, and the demand for high-speed and high-resolution microscopic imaging technology is growing. Especially in the dynamic imaging of biological tissues and cells, traditional optical microscopy techniques often cannot provide imaging speeds and resolutions that meet the needs of modern biomedical research due to inherent limitations of the system. Fluorescence microscopy has become an indispensable tool in biomedical research because of its ability to provide submicron spatial resolution, molecular specificity, and high contrast. These characteristics enable researchers to directly observe and analyze biological structures and functions. However, the inherent blur and noise in the imaging process, especially Poisson noise, often degrade the quality of fluorescence data, resulting in imperfect estimates of the sample.

[0003] To address these issues, deconvolution techniques are widely used in fluorescence microscopy to partially reverse these degradations and improve image resolution and contrast. The Richardson-Lucy deconvolution (RLD) algorithm is an effective method that can use the known point spread function (PSF) and damaged image data to gradually remove the blur of sample density estimation through iterative processing. Each iteration aims to sharpen the image and make it closer to the representation of the real sample.

[0004] Although the RLD algorithm performs well in improving image quality, its computationally intensive nature makes processing speed a limiting factor. In practical applications, especially when large volumes or three-dimensional data sets need to be processed, obtaining high-resolution results often requires multiple iterations, each of which requires significant computing resources and time, and is computationally inefficient, making it difficult to meet real-time or near-real-time imaging requirements. While hoping to increase the number of iterations to obtain high-quality images, image overfitting may occur, that is, artifacts are generated in the image. Therefore, in theory, it is not possible to improve image quality by increasing the number of deconvolution iterations, which will greatly limit the imaging speed and quality.

[0005] In summary, there is currently a lack of a fast three-dimensional deconvolution microscopy method that is fast, has high resolution and imaging contrast, and can achieve three-dimensional imaging. Summary of the invention

[0006] In order to overcome the above-mentioned deficiencies of the prior art, the present invention provides a fast three-dimensional deconvolution microscopy imaging method based on mismatched back projection.

[0007] The technical solution to achieve the purpose of the present invention is: a fast three-dimensional deconvolution microscopy imaging method based on mismatched back projection, the specific steps are:

[0008] Step 1: Take a fluorescent 3D image stack;

[0009] Step 2: Calculate the three-dimensional absorption transfer function and point spread function PSF (forward projection f);

[0010] Step 3: Perform a Wiener-Butterworth filter on the mirror image of the point spread function to generate an unmatched back projection;

[0011] Step 4: Perform fast 3D deconvolution iterations on the 3D image stack using backprojection to deblur the image.

[0012] Preferably, in step 1, a fluorescence incoherent microscopy mode is used to capture a fluorescence three-dimensional image stack, a circular fluorescence light source and a circular pupil are used for illumination, the emission wavelength and excitation wavelength of the fluorescence match the wavelength band of the fluorescent sample, and the camera captures the image in the fluorescence incoherent imaging mode.

[0013] Preferably, the three-dimensional absorption transfer function in step 2 is specifically:

[0014]

[0015] Where μ is the spatial frequency relative to the spatial coordinate x, η is the spatial frequency relative to the spatial coordinate z, S and P are the light source aperture function and pupil function respectively, * represents the complex conjugate, δ is the Dirac function, (x, y, z) is the three-dimensional spatial coordinate of the image, λ is the wavelength of the incident light, and μ′ is the translation of the light source aperture function in the spatial frequency domain.

[0016] Preferably, the point spread function is obtained by inverse Fourier transforming the three-dimensional absorption transfer function, and is expressed as: in is the Fourier transform.

[0017] Preferably, the light source aperture function adopts a circular light source S circ (μ):

[0018]

[0019] Where μ = (μ, ξ) is the spatial frequency relative to the spatial coordinates x and y, and ρ s is the radius of the circle-centered illumination pattern, ρ a and ρ b are the inner and outer radii of the ring pattern, respectively;

[0020] The pupil function P(μ) is:

[0021]

[0022] Where μ = (μ, ξ) is the spatial frequency relative to the spatial coordinates x and y, and ρ p is the radius of the circular pupil function.

[0023] Preferably, in step 3, the specific method of performing Wiener-Butterworth filtering on the mirror image of the point spread function to generate the mismatched back projection is:

[0024] Step 3.1: Flip the point spread function to get the matching backprojector PSF flip , for the back-projector PSF flip Perform fast Fourier transform to obtain the optical transfer function OTF flip ;

[0025] Step 3.2: Optical Transfer Function OTF flip Normalized to get OTF flip_norm , for OTF flip_norm Perform Wiener filtering to obtain OTF Wiener , and calculate the Wiener OTF Wiener The filter cutoff frequency β wienerx ;

[0026] Step 3.3: Set the width k of the Butterworth filter cx , k cy , k cz , adjust the order p of the high-order low-pass filter n , and generate a Butterworth filter;

[0027] Step 3.4: Create the Wiener-Butterworth filter:

[0028]

[0029] Where * represents convolution.

[0030] Step 3.5: Backprojector PSF calculated in step 3.1 flip Using the Wiener-Butterworth filter, an improved mismatched backprojector PSF is generated WB .

[0031] Preferably, the specific formula of Wiener filtering is:

[0032] OTF Wiener =OTF flip_norm / (|OTF flip_norm | 2 +α)

[0033] In the formula, α is taken as 0.05.

[0034] Preferably, the generated Butterworth filter is specifically:

[0035]

[0036] In the formula, p n is the order of the Butterworth high-order filter, w represents the radial frequency in the frequency domain, β represents the OTF value at the limiting resolution;

[0037] Preferably, the specific method of performing fast three-dimensional deconvolution iteration on the three-dimensional image stack by back projection to achieve image deblurring is:

[0038] Step 4.1: preprocessing the three-dimensional image stack, stacking the three-dimensional image stacks at different depths into a three-dimensional array;

[0039] Step 4.2: Convert the entire 3D image dataset into a fluorescent 3D stack object 3d , 3D absorption transfer function, point spread function PSF, mismatched backprojector PSF WB One-time transfer to GPU memory;

[0040] Step 4.3: Run the deconvolution iteration on the GPU, specifically using the Richardson-Lucy deconvolution process, specifically:

[0041] Set the pre-iteration E k The initial value is Object 3d , in each iteration, calculate E k Convolution with PSF Object 3d and Quotient (corrected value) Correction value and the back-projected PSF WB Convolution Combine it with the current estimated value E of the image k Multiply them together to get the updated estimate E of the image k+1 , and used as the initial value E for the next iteration k .

[0042] Preferably, the convolution process is performed by FFT3 and FFTshift3, first calculating E k And PSF data HI and Hpsf in the frequency domain:

[0043] HI = fftshift3(fft3(E k ))

[0044] Hpsf=fftshift3(fft3(ifftshift3(PSF)))

[0045] Recalculate

[0046] Compared with the prior art, the present invention has the following significant advantages: (1) The present invention uses circular partially coherent illumination, and the spatial resolution obtained is close to the diffraction limit under incoherent illumination. (2) The present invention uses a deconvolution method based on mismatched back projection to deblur and suppress noise in images under low light efficiency and high-speed exposure, thereby obtaining high-quality images and greatly improving the speed of the deconvolution algorithm. (3) The present invention uses CUDA parallel processing technology to further optimize RL deconvolution, and utilizes the floating-point computing power of the GPU to increase the speed by another 30-40 times, thereby achieving efficient and fast calculation of large-scale data. (4) The fast three-dimensional deconvolution microscopy imaging method based on mismatched back projection provided by the present invention can realize real-time three-dimensional deconvolution, has broad application prospects, and plays an important role in drug development, clinical diagnosis, etc.

[0047] The present invention is further described in detail below in conjunction with the accompanying drawings and specific embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 The present invention is a flow chart of the fast three-dimensional deconvolution microscopy imaging method based on mismatched back projection.

[0049] Figure 2 It is a USAF resolution target (512x512) and a three-dimensional object (512x512x64) composed of the USAF resolution target tilted in three-dimensional space.

[0050] Figure 3 PSF and traditional backprojection and mismatched backprojection PSF WB Display in real and Fourier space.

[0051] Figure 4 Experimental results of deconvolution of tilted USAF resolution targets under traditional (matched) backprojection and improved unmatched backprojection

[0052] Figure 5 To capture the tilted USAF resolution target (red dashed line), the effects of using traditional backprojection with a high number of iterations and using mismatched backprojection with a low number of iterations are shown.

[0053] Figure 6 These are the deconvolution experimental results of Sphaerotheca sphaeroides under traditional (matched) back-projection and improved unmatched back-projection.

[0054] Figure 7 Statistical effect of accelerating deconvolution iterations using CUDA parallel processing technology. DETAILED DESCRIPTION

[0055] In order to facilitate the understanding and implementation of the present invention, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the implementation examples described herein are only used to illustrate and explain the present invention and are not used to limit the present invention.

[0056] Please see Figure 1 The present invention provides a fast three-dimensional deconvolution microscopy imaging method based on mismatched back projection, comprising the following contents:

[0057] After characterizing the microscopic imaging system, the blurring process of the image can be described and an iterative deconvolution algorithm can be used to reverse this degradation from a noisy, blurred image and estimate the original image. The key components of this method are the forward projector (point spread function PSF), which describes the mapping from the desired image of the object to the noisy, blurred image measured by the microscope, and the back projector, which maps the measured image back to the desired object image. Based on the Richardson-Lucy deconvolution, this method uses the Wiener-Butterworth filter to generate an unmatched back projector to restore the image, which greatly reduces the number of iterations required to restore to the image resolution limit while obtaining a high-resolution image, and reduces image artifacts caused by the increased number of deconvolutions.

[0058] For ease of implementation, the present invention deblurs and suppresses noise for images under low light efficiency and high-speed exposure based on RL deconvolution with mismatched backprojection: The core idea is based on maximum likelihood estimation, aiming to find an image that best explains the observed data. This method assumes that the degradation of the image is linear and spatially invariant, that is, the observed image is the convolution of the ideal image with the point spread function (PSF), plus Poisson noise. The RL algorithm gradually approximates the original unblurred image through an iterative process.

[0059] After choosing an initial estimate, which can usually be the observed image itself, iterations are performed to recover the image. Each iteration consists of two main steps: computing the predicted image: convolving the current estimated image with the forward projector PSF to obtain the predicted blurred image; computing the ratio and updating the image: dividing the observed image by the predicted image to obtain a ratio image; then convolving this ratio image with the flip (or mirror) of the back projector PSF, and using this result to multiply the current image estimate element-wise to update the image estimate.

[0060] In a specific implementation, after characterizing the imaging system, a Richardson-Lucy deconvolution iteration formula is obtained, and the optical transfer function and point spread function of the system are calculated to generate an improved mismatched reverse projector. The steps include:

[0061] Step 1: Take a fluorescent 3D image stack 3d , characterize the microscopic imaging system; use a circular fluorescent light source and a circular pupil, select the emission wavelength and excitation wavelength of the fluorescence excitation module to match the band of the fluorescent sample, and the camera captures images in the fluorescence incoherent imaging mode.

[0062] Step 2: Calculate the three-dimensional absorption transfer function and point spread function PSF (forward projector f) of the system to characterize the image blur caused by the band-limited microscope;

[0063] Among them, the three-dimensional absorption transfer function is specifically:

[0064]

[0065] Wherein, μ = (μ, ξ) is the spatial frequency relative to the spatial coordinates x and y, μ′ = (μ′, ξ′) is the translation of S(μ) and P(μ) in the spatial frequency domain, η is the spatial frequency relative to the spatial coordinate z, S(μ) and P(μ) are the light source aperture function and pupil function, respectively, * represents the complex conjugate, δ is the Dirac function, (x, y, z) is the three-dimensional spatial coordinates of the image, and λ is the wavelength of the incident light.

[0066] Among them, the light source aperture function adopts a circular light source S circ (μ):

[0067]

[0068] Where μ = (μ, ξ) is the spatial frequency relative to the spatial coordinates x and y, and ρ s is the radius of the circle-centered illumination pattern, ρ a and ρ b are the inner and outer radii of the ring pattern respectively.

[0069] The pupil function P(μ) is:

[0070]

[0071] Where μ = (μ, ξ) is the spatial frequency relative to the spatial coordinates x and y, and ρ p is the radius of the circular pupil function.

[0072] Step 3: Mirror PSF of the Point Spread Function PSF flip (matched back projection) performs Wiener-Butterworth filtering to generate an unmatched back projection b. The specific method is:

[0073] Step 3.1: Flip the point spread function PSF (forward projection f) to get the matching backprojector PSF flip , and perform fast Fourier transform to obtain the optical transfer function OTFflip , for the following filtering operation;

[0074] Step 3.2: Normalize the optical transfer function and calculate the normalized optical transfer function OTF. flip_norm Perform Wiener filtering, the specific formula is:

[0075] OTF Wiener =OTF flip_norm / (|OTF flip_norm | 2 +α)

[0076] In the formula, α is taken as 0.05, and the cutoff frequency β of the Wiener filter is calculated wienerx , used for subsequent parameter setting of Butterworth filter;

[0077] Step 3.3: Set the width k of the Butterworth filter cx , k cy , k cz , adjusting the order n of the high-order low-pass filter to 10 to adapt to the specific imaging system characteristics and generate a Butterworth filter;

[0078]

[0079] Step 3.4: Create the Wiener-Butterworth filter:

[0080]

[0081] Where Winer is the Wiener filter, and ee is obtained from the cutoff frequency of the Wiener filter: w represents the radial frequency in the frequency domain, which is calculated by the distance between the current pixel and the center of the image and the width of the Butterworth filter in the corresponding direction (k cx , k cx , k cx ) ratio, p n is the order of the Butterworth high-order filter;

[0082] Step 3.5: PSF calculated in step 3.1 flip Generate an improved mismatched backprojector PSF using the Wiener-Butterworth filter WB .

[0083] Step 4: Perform fast 3D deconvolution iterations on the 3D image stack using backprojection to achieve image deblurring.

[0084] Set the pre-iteration E k The initial value is Object 3d , in each iteration, calculate E k Convolution with PSF Object 3d and Quotient (corrected value) Correction value and the back-projected PSF WB Convolution Combine it with the current estimated value E of the image k Multiply them together to get the updated estimate E of the image k+1 , and used as the initial value E for the next iteration k .

[0085] That is, the specific formula for Richardson-Lucy deconvolution processing is:

[0086]

[0087] Where E k is the current estimate of the desired object image, E k+1 For future estimated value, Object 3d is the measured image, PSF is the forward projection, PSF WB is the back projection, Represents convolution.

[0088] Replacing PSF in Richardson-Lucy deconvolution flip PSF WB After that, the image is reconstructed by deconvolution iteration to improve the iterative convergence speed while ensuring the image quality.

[0089] The present invention is based on the unmatched back-projector and CUDA parallel processing technology, which can realize 3D image reconstruction with high precision and high efficiency, and greatly improve the imaging speed while restoring the image to the theoretical resolution limit, and reduce the artifacts caused by the increase in the number of deconvolution times, and can realize real-time imaging to meet the urgent needs of modern biomedical imaging.

[0090] The present invention performs CUDA acceleration on Richardson-Lucy deconvolution and optimizes serial programs using CUDA parallel processing technology, including: configuring and optimizing the number of cores and thread layout for algorithm execution in a CUDA environment to maximize the parallel processing capability of the GPU; loading the entire three-dimensional data set into the GPU memory to reduce data transmission between the CPU and the GPU, thereby reducing latency and increasing data processing speed; executing the entire iterative deconvolution process on the GPU, including deconvolution operations performed by applying an unmatched backprojector, which utilizes the parallel processing capability of CUDA to accelerate each iterative calculation; and performing FFTSh on the GPU. Ift3d is optimized by calculating the spatial coordinates x, y, z corresponding to each CUDA thread in the three-dimensional data block, and calculating its new spatial coordinates x1, y1, z1 after FFTShift3d. The optimized CUDA kernel function is used to exchange the data of each thread once, without the need to perform multiple operations of one-dimensional FFTShift1d on the Z axis after layer-by-layer FFTShift2d. The CUDA fast Fourier transform library (CUFFT) is used to perform FFT3 forward and inverse transformations. These transformations are the core parts of calculating convolution and deconvolution in the deconvolution process. Compared with the CPU, the FFT3 forward and inverse transformations can be accelerated by 80 times.

[0091] Example

[0092] In order to verify the effectiveness of the present invention, the performance of the imaging method was verified through experiments.

[0093] In order to quantitatively verify the performance of the RL deconvolution method based on the mismatched backprojector, experiments were first conducted on the USAF resolution board. Figure 2 The resolution of the target USAF is tilted in three-dimensional space to form a three-dimensional object.

[0094] Figure 3 Figure a shows the forward projection f in the microscope system. Figure a1 is a transverse slice of the PSF in real space, Figure a2 is its display in the frequency domain, and the right side is the display in the axial slice. Figures b and c are the displays of the Wiener filter and Butterworth filter in real space and frequency domain, where b3 and c3 are the effects of two filtering on the forward projection f (see Figure a1), and Figures d and e are different back projections (bp), including the traditional back projection (PSF transpose) and WB back projection commonly used in RLD. The first two rows are shown in Figure a, and the last row shows the product of the forward and back projections in the frequency domain.

[0095] Figure 4In the middle, a1-a3 are the images of the cross section of the three-dimensional object z=31-33 after systematic blurring. The right side is the comparison of the deconvolution results using the traditional method and the improved mismatched back projector. The improved mismatched back projector is used for 3D RL deconvolution to restore the three-dimensional object for visualization, and the number of iterations can be reduced by 10 times while obtaining the same quality image.

[0096] like Figure 5 As shown in the figure, a segment of the USAF resolution target under a specific XY section is intercepted (the red dashed line part). Figure 5 a, and compare the effects of traditional (Tra) and improved (WB) back-projectors in RL deconvolution, Figure 5 b1-b2 are the comparisons of the traditional method and the improved mismatched backprojector after low and high deconvolution iterations, respectively, comparing the imaging resolution and speed of the two methods. The three-dimensional fluorescence stack (801x801x41) of the spherical algae was tested. Figure 6 As shown, Figure 6 a1-a3 are the cross sections of the original images taken at z=21-23, and the right side shows the processing results using the traditional method and the improved mismatched back projector. When restoring the image to the resolution limit, the number of iterations required can still be reduced by more than 10 times.

Claims

1. A fast three-dimensional deconvolution microscopy imaging method based on mismatched backprojection, characterized in that: The specific steps are: Step 1: Take a fluorescent 3D image stack; Step 2: Calculate the three-dimensional absorption transfer function and point spread function PSF. The three-dimensional absorption transfer function is specifically: Where, μ is the spatial frequency relative to the spatial coordinate x, η is the spatial frequency relative to the spatial coordinate z, S and P are the light source aperture function and pupil function respectively, * represents the complex conjugate, (x, y, z) is the three-dimensional spatial coordinate of the image, λ is the wavelength of the incident light, and μ′ is the translation of the light source aperture function in the spatial frequency domain; The light source aperture function uses a circular light source S circ (μ): Where μ = (μ, ξ) is the spatial frequency relative to the spatial coordinates x and y, and ρ s is the radius of the circle-centered illumination pattern, ρ a and ρ b are the inner and outer radii of the ring pattern, respectively; The pupil function P(μ) is: Where μ = (μ, ξ) is the spatial frequency relative to the spatial coordinates x and y, and ρ p is the radius of the circular pupil function; Step 3: Perform a Wiener-Butterworth filter on the mirror image of the point spread function to generate an unmatched back projection. The specific method is: Step 3.1: Flip the point spread function to get the matching backprojector PSF flip , for the back-projector PSF flip Perform fast Fourier transform to obtain the optical transfer function OTF flip ; Step 3.2: Optical Transfer Function OTF flip Normalized to get OTF flip_norm , for OTF flip_norm Perform Wiener filtering to obtain OTF Wiener , and calculate the Wiener OTF Wiener The filter cutoff frequency β wienerx ; Step 3.3: Set the width k of the Butterworth filter cx , k cy , k cz , adjust the order p of the high-order low-pass filter n , and generate a Butterworth filter; Step 3.4: Create the Wiener-Butterworth filter: In the formula, represents convolution, p n is the order of the Butterworth high-order filter, w represents the radial frequency in the frequency domain, Winer is the Wiener filter, β represents the OTF value at the limit resolution; p n is the order of the Butterworth high-order filter, w represents the radial frequency in the frequency domain; Step 3.5: Backprojector PSF calculated in step 3.1 flip Using the Wiener-Butterworth filter, an improved mismatched backprojector PSF is generated WB ; Step 4: Perform fast 3D deconvolution iterations on the 3D image stack using backprojection to deblur the image as follows: Step 4.1: preprocessing the three-dimensional image stack, stacking the three-dimensional image stacks at different depths into a three-dimensional array; Step 4.2: Convert the entire 3D image dataset into a fluorescent 3D stack object 3d , 3D absorption transfer function, point spread function PSF, mismatched backprojector PSF WB One-time transfer to GPU memory; Step 4.3: Run the deconvolution iteration on the GPU, specifically using the Richardson-Lucy deconvolution process, specifically: Set the pre-iteration E k The initial value is Object 3d , in each iteration, calculate E k Convolution with PSF Object 3d and Business and the back-projected PSF WB Convolution Combine it with the current estimated value E of the image k Multiply them together to get the updated estimate E of the image k+1 , and used as the initial value E for the next iteration k .

2. The fast three-dimensional deconvolution microscopy imaging method based on mismatched back projection according to claim 1, characterized in that: In step 1, a fluorescence three-dimensional image stack is captured using a fluorescence incoherent microscopy imaging mode. A circular fluorescence light source and a circular pupil are used for illumination. The emission wavelength and excitation wavelength of the fluorescence match the wavelength band of the fluorescent sample. The camera captures the image in the fluorescence incoherent imaging mode.

3. The fast three-dimensional deconvolution microscopy imaging method based on mismatched back projection according to claim 1, characterized in that: The point spread function is obtained by inverse Fourier transform of the three-dimensional absorption transfer function and is expressed as: in is the Fourier transform.

4. The fast three-dimensional deconvolution microscopy imaging method based on mismatched back projection according to claim 1, characterized in that: The specific formula of Wiener filtering is: OTF Wiener =OTF flip_norm / (|OTF flip_norm | 2 +α) In the formula, α is taken as 0.

05.

5. The fast three-dimensional deconvolution microscopy imaging method based on mismatched back projection according to claim 1, characterized in that: The generated Butterworth filter is specifically: In the formula, p n is the order of the Butterworth high-order filter, w represents the radial frequency in the frequency domain, and β represents the OTF value at the limiting resolution.

6. The fast three-dimensional deconvolution microscopy imaging method based on mismatched back projection according to claim 1, characterized in that: The convolution process is performed through fft3 and fftshift3. First, E is calculated k And PSF data HI and Hpsf in the frequency domain: HI=fftshift3(fft3(E k )) Hpsf=fftshift3(fft3(ifftshift3(PSF))) Recalculate

Images