A fluorescence image deconvolution method based on multi-scale basis

Abstract

The present invention discloses a fluorescence image deconvolution method based on a multi-scale basis. Starting from the perspective of multi-resolution analysis, the present invention first designs a regularization term based on the characteristics of biological fluorescence images, which have distinct geometric shapes and rich information, and constructs a corresponding deconvolution algorithm, which has better noise reduction and deconvolution capabilities. Furthermore, for time series images, the present invention uses a deconvolution method that first deconvolutes the time domain, then the spatial domain, and then the time domain. Instead of pursuing a one-time noise reduction, the method performs noise reduction and deblurring in steps based on the characteristics of different wavelet bases, achieving both strong noise reduction and better resolution preservation. Furthermore, the present invention utilizes three-dimensional wavelets to fully utilize the three-dimensional continuous information of the sample sequence, combining the advantages and disadvantages of three-dimensional separable transforms and non-separable transforms to design a more optimized noise reduction process.

Description

Technical Field

[0001] The present invention relates to fluorescence microscopy imaging technology, and in particular to a fluorescence image deconvolution method based on a multi-scale basis. Background Art

[0002] In order to ensure that cells can survive for a long time during the imaging process, lower light intensity illumination is usually used in fluorescence microscopy to reduce phototoxicity to cells. In order to capture the rapid activities of cells, shorter exposure times are usually used in fluorescence microscopy. In both cases, the quality of fluorescence microscopy images will eventually deteriorate severely, with a lot of noise, making it difficult to achieve good observation and judgment. Deconvolution post-processing algorithm is an effective technical means to solve this problem. Based on prior knowledge of the imaged biological samples, the deconvolution algorithm can recover the image information that is submerged by noise through computational means. Reference Image reconstruction for structured-illuminationmicroscopy with low signal level [1] A deconvolution algorithm based on global variation regularization (TV) deconvolution is proposed to recover the noisy image information. [2] A deconvolution algorithm based on Hessian deconvolution is proposed to restore noisy image information.

[0003] In traditional deconvolution algorithms, spatial continuity is often exploited as a priori information, so the first-order or second-order derivatives of the image are often minimized as regularization terms. However, as verified in natural image processing, such schemes often only achieve limited optimization results, and using global changes as regularization often inevitably leads to a decrease in resolution.

[0004] Current post-processing deconvolution-based methods for recovering low signal-to-noise ratio data are often based on relatively simple spatial regularization, specifically the spatial continuity of the imaged biological sample, mathematically represented as the derivative of the spatial coordinates. Consequently, the optimization effects achieved are relatively limited. Furthermore, for long-term time series, while existing optimization algorithms utilize changes in the temporal coordinates (i.e., partial derivatives) for optimization, they are essentially still optimizations of two-dimensional slices and fail to truly utilize three-dimensional structural information. Summary of the Invention

[0005] In order to overcome the shortcomings of the above-mentioned existing technologies, the present invention proposes a fluorescence image deconvolution method based on a multi-scale basis, establishes an optimization method for fluorescence images using a multi-scale basis, and utilizes the powerful signal-noise discrimination ability of the multi-scale basis to achieve optimization of fluorescence images with various signal-to-noise ratios, thereby improving and enhancing the existing fluorescence image deconvolution method.

[0006] For a two-dimensional fluorescence microscopic image containing only spatial plane information, the fluorescence image deconvolution method based on a multi-scale basis of the present invention comprises the following steps:

[0007] 1) Obtain original two-dimensional fluorescence images:

[0008] An original two-dimensional fluorescence image is acquired through an optical microscope, and then the original two-dimensional fluorescence image is normalized to obtain a normalized original two-dimensional fluorescence image;

[0009] 2) Deconvolution of the original two-dimensional fluorescence image based on the joint sparse regularization of piecewise linear frame wavelet and curvelet:

[0010] i. Calculating a point spread function for deconvolution using a corresponding theoretical fitting formula based on the wavelength of the illumination light used in the imaging process and the pixel size and numerical aperture of the sensor; or constructing a point spread function for deconvolution using a point spread function actually measured during the imaging process; then normalizing the point spread function to obtain a normalized point spread function;

[0011] ii. Establish a two-dimensional piecewise linear frame wavelet and a two-dimensional curvelet respectively;

[0012] iii. Construct the deconvolution model as follows:

[0013]

[0014] Among them, the first term is the fidelity term, which ensures the fidelity of the deconvolution algorithm, b is the normalized original two-dimensional fluorescence image, f is the image after deconvolution iteration, A is the matrix form of the point spread function, the operator W represents the two-dimensional piecewise linear frame wavelet forward transform, the operator C represents the two-dimensional curvelet forward transform, λ1 and λ2 are two regularization parameters, ||1 and ||2 represent the first norm and the second norm, respectively;

[0015] iv. Solving the above deconvolution model based on the fast soft threshold filtering iteration (FISTA) method to obtain the deconvolution iterative image;

[0016] 3) Airspace optimization:

[0017] For the image after the deconvolution iteration, Richardson-Lucy deconvolution iteration optimization is performed to obtain a deconvolution optimized image.

[0018] In step 1), any form of optical microscope can be used to acquire the original two-dimensional fluorescence image; due to the optimization of the post-processing of the deconvolution algorithm, the quality requirements for the two-dimensional fluorescence microscopic image acquired by the front end can be relaxed when acquiring the image, allowing original images of poor quality.

[0019] In step 2) i., based on the wavelength of the illumination light used in the imaging process and the pixel size and numerical aperture of the sensor, the point spread function used for deconvolution is calculated using the corresponding theoretical fitting formula, using the Bessel fitting formula or the Gaussian formula, etc.

[0020] In step 2) ii., a two-dimensional piecewise linear frame wavelet is established, comprising the following steps:

[0021] a) The coefficients of the frequency domain filter for constructing the piecewise linear frame wavelet are [3]:

[0022]

[0023] Among them, u0 is the coefficient of the low-pass filter, u1 is the coefficient of the first high-pass filter, and u2 is the coefficient of the second high-pass filter;

[0024] b) Performing an inverse Fourier transform on the frequency domain filter of the piecewise linear frame wavelet to the spatial domain to obtain three functions of the one-dimensional piecewise linear frame wavelet, which are expressed as φ(x) or φ(y), ψ1(x) or ψ1(y) and ψ2(x) or ψ2(y) for the two directions x and y of the two-dimensional image respectively;

[0025] c) Extend the one-dimensional piecewise linear frame wavelet to two dimensions through tensor product:

[0026]

[0027]

[0028] Among them, Φ(x,y) is the scaling function, and {Ψ(x,y)} is the wavelet function set, thereby obtaining a two-dimensional spatial domain piecewise linear frame wavelet.

[0029] In step 2) ii., establishing a two-dimensional curve wave includes the following steps:

[0030] a) Construct anisotropic wedge-shaped region U in the frequency domain j [4]:

[0031]

[0032] Where r is the orientation length in polar coordinates, θ is the angle, W is the radial window, and V is the angular window. Indicates the rounding operation, j = 0, 1, 2... is the scale parameter;

[0033] b) For the wedge-shaped region U in the frequency domain j Perform an inverse Fourier transform on (r,θ) to obtain a two-dimensional curve wave.

[0034] In step 2) iii., the two-dimensional piecewise linear frame wavelet forward transform refers to taking the inner product of the two-dimensional image and the two-dimensional piecewise linear frame wavelet to obtain the decomposition coefficients of the two-dimensional piecewise linear frame wavelet. The two-dimensional curvelet forward transform refers to decomposing the two-dimensional image using the wrapping method published in the literature [4] to obtain the two-dimensional curvelet decomposition coefficients.

[0035] In step 2) iv., solving the deconvolution model by the FISTA method to obtain the deconvolution iterative image includes the following steps:

[0036] (1) Calculate the Lipschitz constant of the normalized original two-dimensional fluorescence image;

[0037] (2) Selecting the regularization parameters λ1 and λ2 of the deconvolution model based on the quality of the original two-dimensional fluorescence image; for original two-dimensional fluorescence images with a peak signal-to-noise ratio (PSNR) of 10 dB or above, the values ​​of λ1 and λ2 are selected within the range of 0.0001 to 0.01; for original two-dimensional fluorescence images with a peak signal-to-noise ratio (PSNR) below 10 dB, the values ​​of λ1 and λ2 are selected within the range of 0.01 to 1;

[0038] (3) Based on the fast soft threshold filtering, the deconvolution model is solved iteratively. First, the parameters are initialized, and then the following iterative strategy is used until the loss function convergence:

[0039]

[0040] Among them, f k-1 、f k and f k+1 are the images obtained by the k-1th, kth and k+1th iterations respectively, f k ' and f k ” is the intermediate image in the iterative process, t k is the step size used in the kth iteration, is the gradient operator, L is the normalized original 2D fluorescence image’s Puschitz constant, W and W T Represent the two-dimensional piecewise linear frame wavelet forward transform and inverse transform, C and C respectively T Represent the two-dimensional curve wave forward transform and inverse transform, T λ is the soft threshold filter operator, defined as follows:

[0041] T λ (M)=(|M|-λ)sgn(M)

[0042] Among them, M is the matrix filtered by soft threshold, λ is the soft threshold filter threshold, sgn is the sign function, and the loss function After convergence, the image after deconvolution iteration is obtained;

[0043] (4) Determine whether the image after deconvolution iteration obtained in step (3) meets the actual use requirements. If it does not meet the actual use requirements, adjust the values ​​of λ1 and λ2 accordingly and enter step (5). If it meets the requirements, end and obtain the image after deconvolution iteration;

[0044] (5) After adjusting the values ​​of λ1 and λ2, repeat steps (3) and (4) until the image after deconvolution iteration meets the actual usage requirements, and obtain the image after deconvolution iteration.

[0045] In step (4), if the image after the deconvolution iteration obtained in step (3) does not meet the actual use requirements, the sum of λ1 and λ2 is adjusted based on the current λ1 and λ2 values ​​through the following two steps: a) if the signal-to-noise ratio of the image after the deconvolution iteration obtained in step (3) is lower than the actual use requirements, the λ1 and λ2 values ​​are increased at the same time by 10 to 30%. Conversely, if the resolution of the image after the deconvolution iteration is lower than the actual use requirements, the λ1 and λ2 values ​​are reduced at the same time by 10 to 30%. b) if the continuity of the image after the deconvolution iteration obtained in step (3) does not meet the actual use requirements, the ratio of λ2 to λ1 is increased by 10 to 30%. Conversely, if the image is overly smooth, the ratio of λ2 to λ1 is reduced by 10 to 30%.

[0046] In step 3), the number of iterations is 1 to 5.

[0047] For a two-dimensional fluorescence microscopic image and a three-dimensional fluorescence microscopic image containing a time series, the fluorescence image deconvolution method based on a multi-scale basis of the present invention comprises the following steps:

[0048] 1) Obtain original 3D sequence images:

[0049] Acquiring original two-dimensional fluorescence image time series or original three-dimensional fluorescence microscopic images through an optical microscope, the original two-dimensional fluorescence image time series or the original three-dimensional fluorescence microscopic images are collectively referred to as original three-dimensional sequence images; performing normalization processing on the original three-dimensional sequence images to obtain normalized original three-dimensional sequence images;

[0050] 2) Pre-denoising based on 3D piecewise linear frame wavelet transform:

[0051] A three-dimensional piecewise linear frame wavelet is established, and the original three-dimensional sequence image is subjected to a three-dimensional piecewise linear frame wavelet forward transform to obtain the decomposition coefficients of the three-dimensional piecewise linear frame wavelet; the decomposition coefficients are subjected to a soft threshold filter and then an inverse three-dimensional piecewise linear frame wavelet transform to obtain the pre-denoised three-dimensional sequence image;

[0052] 3) Treat each of the pre-denoised 3D sequence images as the original 2D fluorescence image and perform deconvolution on them one by one based on the joint sparse regularization of piecewise linear frame wavelet and curvelet:

[0053] i. normalizing the original two-dimensional fluorescence image to obtain a normalized original two-dimensional fluorescence image;

[0054] ii. Based on the wavelength of the illumination light used in the imaging process and the pixel size and numerical aperture of the sensor, a point spread function (PSF) for deconvolution is calculated using a corresponding theoretical formula; alternatively, a PSF actually measured during the imaging process is used to construct the PSF for deconvolution; the PSF is then normalized to obtain a normalized PSF;

[0055] iii. Establish a two-dimensional piecewise linear frame wavelet and a two-dimensional curvelet respectively;

[0056] iv. Construct the deconvolution model as follows:

[0057]

[0058] Among them, the first term is the fidelity term, which ensures the fidelity of the deconvolution algorithm, b is the normalized original two-dimensional fluorescence image, f is the image after deconvolution iteration, A is the matrix form of the point spread function, the operator W represents the two-dimensional piecewise linear frame wavelet forward transform, the operator C represents the two-dimensional curvelet forward transform, λ1 and λ2 are two regularization parameters, ||1 and ||2 represent the first norm and the second norm, respectively;

[0059] v. Solving the above deconvolution model based on the fast soft threshold filtering iteration (FISTA) method to obtain the deconvolution iterative image;

[0060] vi. Recombining each deconvolution iterative image in the original order to obtain a three-dimensional sequence image after deconvolution iterative;

[0061] 4) Airspace optimization:

[0062] For the three-dimensional sequence images after deconvolution iteration, Richardson-Lucy deconvolution iteration optimization is performed to obtain deconvolution optimized three-dimensional sequence images;

[0063] 5) Three-dimensional optimization based on three-dimensional dual-tree complex wavelet:

[0064] The deconvolution optimized three-dimensional sequence image is subjected to a three-dimensional dual-tree complex wavelet forward transform to obtain its three-dimensional dual-tree complex wavelet decomposition coefficients; the decomposition coefficients are subjected to a soft threshold filter and then a three-dimensional dual-tree complex wavelet inverse transform to obtain a continuity optimized three-dimensional sequence deconvolution image.

[0065] In step v. of step 3), solving the deconvolution model by the FISTA method to obtain the image after deconvolution iteration includes the following steps:

[0066] (1) Calculate the Lipschitz constant of the normalized original two-dimensional fluorescence image;

[0067] (2) Selecting the regularization parameters λ1 and λ2 of the deconvolution model based on the quality of the original two-dimensional fluorescence image; for original two-dimensional fluorescence images with a peak signal-to-noise ratio (PSNR) of 10 dB or above, the values ​​of λ1 and λ2 are selected within the range of 0.0001 to 0.01; for original two-dimensional fluorescence images with a peak signal-to-noise ratio (PSNR) below 10 dB, the values ​​of λ1 and λ2 are selected within the range of 0.01 to 1;

[0068] (3) Based on the fast soft threshold filtering, the deconvolution model is solved iteratively. First, the parameters are initialized, and then the following iterative strategy is used until the loss function convergence:

[0069]

[0070] Among them, f k-1 、f k and f k+1 are the images obtained by the k-1th, kth and k+1th iterations respectively, f k ' and f k ” is the intermediate image in the iterative process, t k is the step size used in the kth iteration, is the gradient operator, L is the normalized original 2D fluorescence image’s Puschitz constant, W and W T Represent the two-dimensional piecewise linear frame wavelet forward transform and inverse transform, C and C respectively T Represent the two-dimensional curve wave forward transform and inverse transform, T λ is the soft threshold filter operator, defined as follows:

[0071] T λ (M)=(|M|-λ)sgn(M)

[0072] Among them, M is the matrix filtered by soft threshold, λ is the soft threshold filter threshold, sgn is the sign function, and the loss function After convergence, the image after deconvolution iteration is obtained;

[0073] (4) Determine whether the image after deconvolution iteration obtained in step (3) meets the actual use requirements. If it does not meet the actual use requirements, adjust the values ​​of λ1 and λ2 accordingly and enter step (5). If it meets the requirements, end and obtain the image after deconvolution iteration;

[0074] (5) After adjusting the values ​​of λ1 and λ2, repeat steps (3) and (4) until the image after deconvolution iteration meets the actual usage requirements, and obtain the image after deconvolution iteration.

[0075] In step (4), if the image after the deconvolution iteration obtained in step (3) does not meet the actual use requirements, the sum of λ1 and λ2 is adjusted based on the current λ1 and λ2 values ​​through the following two steps: a) if the signal-to-noise ratio of the image after the deconvolution iteration obtained in step (3) is lower than the actual use requirements, the λ1 and λ2 values ​​are increased at the same time by 10 to 30%. Conversely, if the resolution of the image after the deconvolution iteration is lower than the actual use requirements, the λ1 and λ2 values ​​are reduced at the same time by 10 to 30%. b) if the continuity of the image after the deconvolution iteration obtained in step (3) does not meet the actual use requirements, the ratio of λ2 to λ1 is increased by 10 to 30%. Conversely, if the image is overly smooth, the ratio of λ2 to λ1 is reduced by 10 to 30%.

[0076] In step 5), the three-dimensional dual-tree complex wavelet forward transform refers to decomposing the three-dimensional sequence image using the method in the literature [5] to obtain the decomposition coefficients of its three-dimensional dual-tree complex wavelet.

[0077] Advantages of the present invention:

[0078] For the first time, the present invention starts from the perspective of multi-resolution analysis and designs regularization terms based on the piecewise linear frame wavelet and curve wave that are sensitive to image edge and direction information according to the characteristics of biological fluorescence images with clear geometric shapes, and constructs a corresponding deconvolution algorithm with better noise reduction and deblurring capabilities. At the same time, for time series images, the present invention uses a deconvolution method of first in the time domain, then in the spatial domain, and then in the time domain. It does not pursue one-time noise reduction, but performs noise reduction and deblurring in steps according to the characteristics of different wavelet bases, which has stronger noise reduction capabilities and better resolution retention capabilities. At the same time, the three-dimensional continuous information of the sample sequence is fully utilized through the three-dimensional wavelet, and the advantages and disadvantages of the three-dimensional piecewise linear frame wavelet, a three-dimensional separable transform, and the three-dimensional dual-tree complex wavelet, an inseparable transform, are combined to design a better noise reduction process. BRIEF DESCRIPTION OF THE DRAWINGS

[0079] Figure 1 Flowchart of Example 1 of the fluorescence image deconvolution method based on multi-scale basis of the present invention;

[0080] Figure 2This is a flow chart of Example 2 of the fluorescence image deconvolution method based on a multi-scale basis of the present invention;

[0081] Figure 3 Figure 1 is an image of actin and its spectrum obtained according to an embodiment of the fluorescence image deconvolution method based on a multiscale basis of the present invention, wherein (a) is an original two-dimensional fluorescence image of actin filaments captured using a structured light illumination microscope and its spectrum, (b) is an image and its spectrum obtained using traditional Hessian deconvolution, and (c) is an image and its spectrum obtained using the fluorescence image deconvolution method based on a multiscale basis proposed in the present invention;

[0082] Figure 4 The present invention provides a long-term endoplasmic reticulum image obtained by an embodiment of the fluorescence image deconvolution method based on a multi-scale basis, wherein (a) is a three-dimensional original sequence image of the endoplasmic reticulum captured by a wide-field microscope, (b) is a pre-denoised three-dimensional original sequence image obtained by pre-denoising using a three-dimensional piecewise linear frame wavelet, (c) is a deconvolution-optimized three-dimensional sequence image obtained after image-by-image deconvolution and spatial domain optimization, and (d) is a final optimized three-dimensional sequence image obtained after three-dimensional complex double-tree wavelet soft threshold filtering. DETAILED DESCRIPTION

[0083] The present invention will be further described below through specific embodiments in conjunction with the accompanying drawings.

[0084] Example 1

[0085] This embodiment processes a two-dimensional fluorescence microscopic image of actin containing only spatial plane information. The embodiment uses a fluorescence image deconvolution method based on a multi-scale basis, such as Figure 1 As shown, the following steps are included:

[0086] 1) Obtain original two-dimensional fluorescence images:

[0087] The original two-dimensional fluorescence image was obtained by optical microscopy. Figure 3 As shown in (a), the original two-dimensional fluorescence image is then normalized to obtain a normalized original two-dimensional fluorescence image;

[0088] 2) Deconvolution of the original two-dimensional fluorescence image based on the joint sparse regularization of piecewise linear frame wavelet and curvelet:

[0089] i. Based on the wavelength of the illumination light used in the imaging process and the pixel size and numerical aperture of the sensor, a theoretical fitting formula is used to calculate the point spread function used for deconvolution. In this example, the illumination wavelength used is 525 nm, the sensor pixel size is 65 nm, and the numerical aperture is 1.49. The point spread function used for deconvolution is fitted using the Bessel formula;

[0090] ii. Establish a two-dimensional piecewise linear frame wavelet and a two-dimensional curvelet respectively;

[0091] iii. Solve the above deconvolution model based on the fast soft threshold filtering iteration (FISTA) method to obtain the image after deconvolution iteration:

[0092] (1) Calculate the Lipschitz constant of the normalized original two-dimensional fluorescence image;

[0093] (2) The regularization parameters λ1 and λ2 of the deconvolution model were selected based on the quality of the original two-dimensional fluorescence image; the λ1 value was selected to be 0.001 and the λ2 value to be 0.001 based on the signal-to-noise ratio of this image.

[0094] (3) Based on the fast soft threshold filtering iterative solution of the deconvolution model, the parameters are first initialized:

[0095]

[0096] Where f0 is the initial image of the iteration, that is, the original two-dimensional fluorescence image, and t0 is the initial step size in the iteration process. The following iterative strategy is then used until the loss function convergence:

[0097]

[0098] where f k-1 、f k and f k+1 are the images obtained by the k-1th, kth and k+1th iterations respectively, f k ' and f k ” is the intermediate image in the iterative process, t k is the step size used in the kth iteration, is the gradient operator, L is the Puschitz constant of the original two-dimensional fluorescence image, W and W T Represent the two-dimensional piecewise linear frame wavelet forward transform and inverse transform, C and C respectively T Represent the two-dimensional curve wave forward transform and inverse transform, T λ is the soft threshold filter operator, defined as follows:

[0099] T λ (M)=(|M|-λ)sgn(M)

[0100] Where M is the matrix filtered by soft threshold, λ is the soft threshold filter threshold, sgn is the sign function, and the loss function After convergence, the image after deconvolution iteration is obtained;

[0101] (4) Determine whether the image after the deconvolution iteration obtained in step (3) meets the actual use requirements. If the image after the deconvolution iteration obtained in step (3) does not meet the actual use requirements, adjust the sum of λ1 and λ2 through the following two steps based on the current λ1 and λ2 value selection: a) If the signal-to-noise ratio of the image after the deconvolution iteration obtained in step (3) is lower than the actual use requirements, increase the λ1 and λ2 values ​​by 20% at the same time. Conversely, if the resolution of the image after the deconvolution iteration is lower than the actual use requirements, reduce the λ1 and λ2 values ​​by 20% at the same time; b) If the continuity of the image after the deconvolution iteration obtained in step (3) does not meet the actual use requirements, increase the ratio of λ2 to λ1 by 20%. Conversely, if the image is over-smoothed, reduce the ratio of λ2 to λ1 by 20%.

[0102] (5) After adjusting the values ​​of λ1 and λ2, repeat steps (3) and (4) until the image after deconvolution iteration meets the actual use requirements. The deconvolution image is obtained. The λ1 value is selected to be 0.0003, and the λ2 value is selected to be 0.0003.

[0103] 3) Airspace optimization:

[0104] For the image after deconvolution iteration, Richardson-Lucy deconvolution iteration optimization is performed with 2 iterations to obtain the deconvolution optimized image, as shown in Figure 3 (c) shown.

[0105] from Figure 3 It can be seen that the Figure 3 The deconvolution optimized image shown in (c) has a significantly higher signal-to-noise ratio and resolution than Figure 3 (b) The traditional deconvolution image is shown.

[0106] Example 2

[0107] This embodiment processes a two-dimensional fluorescence microscopic image containing a time series. The fluorescence image deconvolution method based on a multi-scale basis in this embodiment is as follows: Figure 2 As shown, the following steps are included:

[0108] 1) Obtain original 3D sequence images:

[0109] The original two-dimensional fluorescence images containing time series are collected by optical microscopy as original three-dimensional sequence images, such as Figure 4 (a)

[0110] 2) Pre-denoising based on 3D piecewise linear frame wavelet transform:

[0111] The original three-dimensional sequence image is subjected to a three-dimensional piecewise linear frame wavelet forward transform to obtain the decomposition coefficients of the three-dimensional piecewise linear frame wavelet; the decomposition coefficients are subjected to a soft threshold filter and then a three-dimensional piecewise linear frame wavelet inverse transform to obtain the three-dimensional sequence image after pre-noise reduction. In this embodiment, the soft threshold filter threshold is 0.1, and the obtained pre-noise reduction image is as follows: Figure 4 (b)

[0112] 3) Treat each of the pre-denoised 3D sequence images as the original 2D fluorescence image and perform deconvolution on them one by one based on the joint sparse regularization of piecewise linear frame wavelet and curvelet:

[0113] i. normalizing the original two-dimensional fluorescence image to obtain a normalized original two-dimensional fluorescence image;

[0114] ii. Based on the wavelength of the illumination light used in the imaging process and the pixel size and numerical aperture of the sensor, a theoretical formula is used to calculate the point spread function used for deconvolution. In this example, the illumination wavelength used is 525 nm, the sensor pixel size is 65 nm, and the numerical aperture is 1.49. The point spread function used for deconvolution is fitted using the Bessel formula.

[0115] iii. Establish a two-dimensional piecewise linear frame wavelet and a spatial curve wavelet respectively;

[0116] iv. Construct the deconvolution model as follows:

[0117]

[0118] Among them, the first term is the fidelity term, which ensures the fidelity of the deconvolution algorithm, b is the normalized original two-dimensional fluorescence image, f is the image after deconvolution iteration, A is the matrix form of the point spread function, the operator W represents the two-dimensional piecewise linear frame wavelet forward transform, the operator C represents the two-dimensional curvelet forward transform, λ1 and λ2 are two regularization parameters, ||1 and ||2 represent the first norm and the second norm, respectively;

[0119] v. Solve the above deconvolution model based on the fast soft threshold filtering iteration (FISTA) method to obtain the image after deconvolution iteration:

[0120] (1) Calculate the Lipschitz constant of the normalized original two-dimensional fluorescence image;

[0121] (2) The regularization parameters λ1 and λ2 of the deconvolution model were selected based on the quality of the original two-dimensional fluorescence image; the regularization parameter of the piecewise linear frame wavelet was selected to be 0.005, and the regularization parameter of the curve wavelet was selected to be 0.005;

[0122] (3) Based on the fast soft threshold filtering iterative solution of the deconvolution model, the parameters are first initialized:

[0123]

[0124] Where f0 is the initial image of the iteration, that is, the original two-dimensional fluorescence image, and t0 is the initial step size in the iteration process.

[0125] Then the following iterative strategy is adopted until the loss function convergence:

[0126]

[0127] where f k-1 、f k and f k+1 are the images obtained by the k-1th, kth and k+1th iterations respectively, f k ' and f k ” is the intermediate image in the iterative process, t k is the step size used in the kth iteration, is the gradient operator, L is the normalized original 2D fluorescence image’s Puschitz constant, W and W T Represent the two-dimensional piecewise linear frame wavelet forward transform and inverse transform, C and C respectively T Represent the two-dimensional curve wave forward transform and inverse transform, T λ is the soft threshold filter operator, defined as follows:

[0128] T λ (M)=(|M|-λ)sgn(M)

[0129] Where M is the matrix filtered by soft threshold, λ is the soft threshold filter threshold, sgn is the sign function, and the loss function After convergence, the image after deconvolution iteration is obtained;

[0130] (4) Determine whether the image after the deconvolution iteration obtained in step (3) meets the actual use requirements. If the image after the deconvolution iteration obtained in step (3) does not meet the actual use requirements, adjust the sum of λ1 and λ2 through the following two steps based on the current λ1 and λ2 value selection: a) If the signal-to-noise ratio of the image after the deconvolution iteration obtained in step (3) is lower than the actual use requirements, increase the λ1 and λ2 values ​​by 20% at the same time. Conversely, if the resolution of the image after the deconvolution iteration is lower than the actual use requirements, reduce the λ1 and λ2 values ​​by 20% at the same time; b) If the continuity of the image after the deconvolution iteration obtained in step (3) does not meet the actual use requirements, increase the ratio of λ2 to λ1 by 20%. Conversely, if the image is over-smoothed, reduce the ratio of λ2 to λ1 by 20%.

[0131] (5) After adjusting the values ​​of λ1 and λ2, repeat steps (3) and (4) until the image after deconvolution iteration meets the actual use requirements, and obtain the image after deconvolution iteration; select the value of λ1 to be 0.005 and the value of λ2 to be 0.002;

[0132] vi. Recombining each deconvolution iterative image in the original order to obtain a three-dimensional sequence image after deconvolution iterative;

[0133] 4) Airspace optimization:

[0134] For the three-dimensional sequence images after deconvolution iteration, Richardson-Lucy deconvolution iteration optimization is performed to obtain deconvolution optimized three-dimensional sequence images. In this embodiment, the number of iterations is 5, and the final deconvolution optimized three-dimensional sequence images are obtained as shown in FIG. Figure 4 (c)

[0135] 5) Time series optimization based on three-dimensional complex dual-tree wavelet:

[0136] The deconvolution optimized three-dimensional sequence image is subjected to a three-dimensional dual-tree complex wavelet forward transform to obtain its three-dimensional dual-tree complex wavelet decomposition coefficients; the decomposition coefficients are subjected to a soft threshold filter and then a three-dimensional dual-tree complex wavelet inverse transform to obtain the final optimized three-dimensional sequence image. In this embodiment, the soft threshold filter threshold is 0.025, and the continuity optimized three-dimensional sequence deconvolution image is obtained as shown in FIG. Figure 4 (d) shown.

[0137] Finally, it should be noted that the purpose of disclosing the embodiments is to facilitate a further understanding of the present invention. However, those skilled in the art will appreciate that various substitutions and modifications are possible without departing from the spirit and scope of the present invention and the appended claims. Therefore, the present invention should not be limited to the contents disclosed in the embodiments; the scope of protection claimed by the present invention shall be determined by the scope defined in the claims.

[0138] References

[0139] [1]Chu K,McMillan PJ,Smith ZJ,et al.Image reconstruction forstructured-illumination microscopy with low signal level[J].Optics Express,2014,22(7):8687-8702.

[0140] [2]Huang X,Fan J,Li L,et al.Fast,long-term,super-resolution imagingwith Hessian structured illumination microscopy[J].Nature Biotechnology,2018,36(5):451-459.

[0141] [3]Daubechies I,Han B,Ron A,et al.Framelets:MRA-based constructionsof wavelet frames[J].Applied and computational harmonic analysis,2003,14(1):1-46.

[0142] [4]Candes E,Demanet L,Donoho D,et al.Fast discrete curvelettransforms[J].multiscale modeling&simulation,2006,5(3):861-899.

[0143] [5]Selesnick I W,Li K Y.Video denoising using 2D and 3D dual-treecomplex wavelet transforms[C] / / Wavelets:Applications in Signal and ImageProcessing X.SPIE,2003,5207:607-618.

Claims

1. A fluorescence image deconvolution method based on a multi-scale basis for a two-dimensional fluorescence microscopy image containing only spatial plane information, characterized in that: The fluorescence image deconvolution method comprises the following steps: 1) Obtain original two-dimensional fluorescence images: An original two-dimensional fluorescence image is acquired through an optical microscope, and then the original two-dimensional fluorescence image is normalized to obtain a normalized original two-dimensional fluorescence image; 2) Deconvolution of the original two-dimensional fluorescence image based on the joint sparse regularization of piecewise linear frame wavelet and curvelet: i. Based on the wavelength of the illumination light used in the imaging process and the pixel size and numerical aperture of the sensor, the point spread function for deconvolution is calculated using the corresponding theoretical fitting formula; or the point spread function actually measured during the imaging process is used to construct the point spread function for deconvolution; the point spread function is then normalized. Get the normalized point spread function; ii. Establish a two-dimensional piecewise linear frame wavelet and a two-dimensional curvelet respectively; iii. Construct the deconvolution model as follows: Among them, the first term is the fidelity term, which ensures the fidelity of the deconvolution algorithm, b is the normalized original two-dimensional fluorescence image, f is the image after deconvolution iteration, A is the matrix form of the point spread function, the operator W represents the two-dimensional piecewise linear frame wavelet forward transform, the operator C represents the two-dimensional curvelet forward transform, λ1 and λ2 are two regularization parameters, ||1 and ||2 represent the first norm and the second norm, respectively; iv. Solving the above deconvolution model based on the fast soft threshold filtering iterative FISTA method to obtain the deconvolution iterative image; 3) Airspace optimization: For the image after the deconvolution iteration, Richardson-Lucy deconvolution iteration optimization is performed to obtain a deconvolution optimized image.

2. The fluorescence image deconvolution method according to claim 1, wherein: In step 2) ii., establishing a two-dimensional piecewise linear frame wavelet comprises the following steps: a) The coefficients of the frequency domain filter for constructing the piecewise linear frame wavelet are: Among them, u0 is the coefficient of the low-pass filter, u1 is the coefficient of the first high-pass filter, and u2 is the coefficient of the second high-pass filter; b) Performing an inverse Fourier transform on the frequency domain filter of the piecewise linear frame wavelet to the spatial domain to obtain three functions of the one-dimensional piecewise linear frame wavelet, which are expressed as φ(x) or φ(y), ψ1(x) or ψ1(y) and ψ2(x) or ψ2(y) for the two directions x and y of the two-dimensional image respectively; c) Extend the one-dimensional piecewise linear frame wavelet to two dimensions through tensor product: Among them, Φ(x,y) is the scaling function, and {Ψ(x,y)} is the wavelet function set, thereby obtaining a two-dimensional spatial domain piecewise linear frame wavelet.

3. The fluorescence image deconvolution method according to claim 1, wherein: In step 2) ii., establishing a two-dimensional curve wave comprises the following steps: a) Construct anisotropic wedge-shaped region U in the frequency domain j : Where r is the orientation length in polar coordinates, θ is the angle, W is the radial window, and V is the angular window. Indicates rounding operation, j = 0, 1, 2... is the scale parameter; b) For the wedge-shaped region U in the frequency domain j Perform an inverse Fourier transform on (r,θ) to obtain a two-dimensional curve wave.

4. The fluorescence image deconvolution method according to claim 1, wherein: In step 2) iv., solving the deconvolution model by the fast FISTA method to obtain the deconvolution iterative image includes the following steps: (1) Calculate the Puschitz constant of the normalized original two-dimensional fluorescence image; (2) Select the regularization parameters λ1 and λ2 of the deconvolution model based on the quality of the original two-dimensional fluorescence image; (3) Based on the fast soft threshold filtering, the deconvolution model is solved iteratively. First, the parameters are initialized, and then the following iterative strategy is used until the loss function convergence: Among them, f k-1 、f k and f k+1 are the images obtained by the k-1th, kth and k+1th iterations respectively, f k ' and f k ” is the intermediate image in the iterative process, t k is the step size used in the kth iteration, is the gradient operator, L is the normalized original 2D fluorescence image’s Puschitz constant, W and W T Represent the two-dimensional piecewise linear frame wavelet forward transform and inverse transform, C and C respectively T Represent the two-dimensional curve wave forward transform and inverse transform, T λ is the soft threshold filter operator, defined as follows: T λ (M)=(|M|-λ)sgn(M) Among them, M is the matrix filtered by soft threshold, λ is the soft threshold filter threshold, sgn is the sign function, and the loss function After convergence, the image after deconvolution iteration is obtained; (4) Determine whether the image after deconvolution iteration obtained in step (3) meets the actual use requirements. If it does not meet the actual use requirements, adjust the values ​​of λ1 and λ2 accordingly and enter step (5). If it meets the requirements, end and obtain the image after deconvolution iteration; (5) After adjusting the values ​​of λ1 and λ2, repeat steps (3) and (4) until the image after deconvolution iteration meets the actual usage requirements, and obtain the image after deconvolution iteration.

5. The fluorescence image deconvolution method according to claim 4, wherein: In step (4), if the image after the deconvolution iteration obtained in step (3) does not meet the actual use requirements, the sum of λ1 and λ2 is adjusted based on the current λ1 and λ2 values ​​through the following two steps: a) if the signal-to-noise ratio of the image after the deconvolution iteration obtained in step (3) is lower than the actual use requirements, the values ​​of λ1 and λ2 are increased at the same time; conversely, if the resolution of the image after the deconvolution iteration is lower than the actual use requirements, the values ​​of λ1 and λ2 are reduced at the same time; b) if the continuity of the image after the deconvolution iteration obtained in step (3) does not meet the actual use requirements, the ratio of λ2 to λ1 is increased; conversely, if the image is overly smooth, the ratio of λ2 to λ1 is reduced.

6. A fluorescence image deconvolution method based on a multi-scale basis for two-dimensional fluorescence microscopic images and three-dimensional fluorescence microscopic images containing time series, characterized in that: The fluorescence image deconvolution method comprises the following steps: 1) Obtain original 3D sequence images: Acquiring original two-dimensional fluorescence image time series or original three-dimensional fluorescence microscopic images through an optical microscope, the original two-dimensional fluorescence image time series or the original three-dimensional fluorescence microscopic images are collectively referred to as original three-dimensional sequence images; performing normalization processing on the original three-dimensional sequence images to obtain normalized original three-dimensional sequence images; 2) Pre-denoising based on 3D piecewise linear frame wavelet transform: A three-dimensional piecewise linear frame wavelet is established, and the original three-dimensional sequence image is subjected to a three-dimensional piecewise linear frame wavelet forward transform to obtain the decomposition coefficients of the three-dimensional piecewise linear frame wavelet; the decomposition coefficients are subjected to a soft threshold filter and then an inverse three-dimensional piecewise linear frame wavelet transform to obtain the pre-denoised three-dimensional sequence image; 3) Treat each of the pre-denoised 3D sequence images as the original 2D fluorescence image and perform deconvolution on them one by one based on the joint sparse regularization of piecewise linear frame wavelet and curvelet: i. normalizing the original two-dimensional fluorescence image to obtain a normalized original two-dimensional fluorescence image; ii. Based on the wavelength of the illumination light used in the imaging process and the pixel size and numerical aperture of the sensor, a point spread function (PSF) for deconvolution is calculated using a corresponding theoretical formula; alternatively, a PSF actually measured during the imaging process is used to construct the PSF for deconvolution; the PSF is then normalized to obtain a normalized PSF; iii. Establish a two-dimensional piecewise linear frame wavelet and a two-dimensional curvelet respectively; iv. Construct the deconvolution model as follows: Among them, the first term is the fidelity term, which ensures the fidelity of the deconvolution algorithm, b is the normalized original two-dimensional fluorescence image, f is the image after deconvolution iteration, A is the matrix form of the point spread function, the operator W represents the two-dimensional piecewise linear frame wavelet forward transform, the operator C represents the two-dimensional curvelet forward transform, λ1 and λ2 are two regularization parameters, ||1 and ||2 represent the first norm and the second norm, respectively; v. Solving the above deconvolution model based on a fast soft threshold filtering iterative method to obtain an image after deconvolution iteration; vi. Recombining each deconvolution iterative image in the original order to obtain a three-dimensional sequence image after deconvolution iterative; 4) Airspace optimization: For the three-dimensional sequence images after deconvolution iteration, Richardson-Lucy deconvolution iteration optimization is performed to obtain deconvolution optimized three-dimensional sequence images; 5) Three-dimensional optimization based on three-dimensional dual-tree complex wavelet: The deconvolution optimized three-dimensional sequence image is subjected to a three-dimensional dual-tree complex wavelet forward transform to obtain its three-dimensional dual-tree complex wavelet decomposition coefficients; the decomposition coefficients are subjected to a soft threshold filter and then a three-dimensional dual-tree complex wavelet inverse transform to obtain a continuity optimized three-dimensional sequence deconvolution image.

7. The fluorescence image deconvolution method according to claim 6, wherein: In step v. of step 3), solving the deconvolution model by the FISTA method to obtain the image after deconvolution iteration includes the following steps: (1) Calculate the Puschitz constant of the normalized original two-dimensional fluorescence image; (2) Select the regularization parameters λ1 and λ2 of the deconvolution model based on the quality of the original two-dimensional fluorescence image; (3) Based on the fast soft threshold filtering, the deconvolution model is solved iteratively. First, the parameters are initialized, and then the following iterative strategy is used until the loss function convergence: Among them, f k-1 、f k and f k+1 are the images obtained by the k-1th, kth and k+1th iterations respectively, f k ' and f k ” is the intermediate image in the iterative process, t k is the step size used in the kth iteration, is the gradient operator, L is the Puschitz constant of the original two-dimensional fluorescence image, W and W T Represent the two-dimensional piecewise linear frame wavelet forward transform and inverse transform, C and C respectively T Represent the two-dimensional curve wave forward transform and inverse transform, T λ is the soft threshold filter operator, defined as follows: T λ (M)=(|M|-λ)sgn(M) Among them, M is the matrix filtered by soft threshold, λ is the soft threshold filter threshold, sgn is the sign function, and the loss function After convergence, the image after deconvolution iteration is obtained; (4) Determine whether the image after deconvolution iteration obtained in step (3) meets the actual use requirements. If it does not meet the actual use requirements, adjust the values ​​of λ1 and λ2 accordingly and enter step (5). If it meets the requirements, end and obtain the image after deconvolution iteration; (5) After adjusting the values ​​of λ1 and λ2, repeat steps (3) and (4) until the image after deconvolution iteration meets the actual usage requirements, and obtain the image after deconvolution iteration.

8. The fluorescence image deconvolution method according to claim 7, wherein: In step (4), if the image after the deconvolution iteration obtained in step (3) does not meet the actual use requirements, the sum of λ1 and λ2 is adjusted based on the current λ1 and λ2 values ​​through the following two steps: a) if the signal-to-noise ratio of the image after the deconvolution iteration obtained in step (3) is lower than the actual use requirements, the values ​​of λ1 and λ2 are increased at the same time; conversely, if the resolution of the image after the deconvolution iteration is lower than the actual use requirements, the values ​​of λ1 and λ2 are reduced at the same time; b) if the continuity of the image after the deconvolution iteration obtained in step (3) does not meet the actual use requirements, the ratio of λ2 to λ1 is increased; conversely, if the image is overly smooth, the ratio of λ2 to λ1 is reduced.

Images