Fast sparse inverse convolution two-dimensional image restoration method

Abstract

The invention discloses a rapid sparse inverse convolution two-dimensional image restoration method. The method comprises the following steps: (1) reading an observation image and a convolution template; (2) carrying out linear transformation on a value range of an observation image pixel, and carrying out normalization processing on a convolution template signal; and (3) constructing a sparse inverse convolution two-dimensional image restoration optimization model, and performing hierarchical and staged rapid sparse inverse convolution restoration. In the first stage, when the L1 norm is minimized, a uniform adjusting parameter lambda is adopted to carry out threshold convergence, and in the second stage, different adjusting parameters are adopted to carry out L1 norm minimization in combination with the weighted L1 norm. According to the method, the L2 norm in the model is minimized through the gradient projection algorithm in combination with the two-dimensional finite discrete fast spatial convolution, then the L1 norm is minimized through the iterative threshold reduction algorithm, and the restored image is obtained through alternate iteration fast sparse solution of the two methods. According to the method provided by the invention, the spatial resolution of the image can be effectively improved, the processing speed is high, and the anti-noise capability is strong.

Description

Technical Field

[0001] This invention belongs to the field of image restoration processing, and more particularly relates to an image deconvolution restoration method. Specifically, it provides a fast sparse deconvolution restoration method for two-dimensional images. The invention is particularly useful for improving the resolution of microscopic optical images. Background Art

[0002] Image restoration is the process of improving the quality of degraded images by removing or mitigating the image degradation that occurs during the digital image acquisition process, thereby achieving a visually improved image. The most typical degradation phenomena are blur and noise. Deconvolution methods can significantly improve image resolution and are therefore widely used in microscopic optical imaging. With the continuous advancement of high-throughput microscopy, large amounts of data can be accumulated in a short period of time, creating an urgent need for fast deconvolution methods to analyze this data. Currently, the principle behind fast deconvolution is that convolution of a signal in the time domain is equivalent to a dot product operation in the frequency domain. Accordingly, deconvolution of a signal in the time domain is equivalent to a dot division operation in the frequency domain. The major drawback of this method is its poor noise immunity. Here, we develop a new fast sparse deconvolution method for 2D image restoration. This method relies on an efficient sparse matrix product operation, avoiding frequency domain analysis of the image. Furthermore, it incorporates prior information about the image. The deconvolution-restored image is closer to the true image and exhibits stronger noise immunity than traditional deconvolution methods. Summary of the Invention

[0003] The present invention aims to provide a fast sparse deconvolution restoration method for two-dimensional images to improve the spatial resolution of the image. This method is characterized by fast processing speed and strong noise immunity. To solve the above technical problems, the present invention adopts the following technical solutions:

[0004] The present invention provides a fast sparse deconvolution two-dimensional image restoration method, comprising the following steps:

[0005] Step 1: Data input. Read the observed image and convolution template.

[0006] The second step is to linearly transform the value range of the observed image pixels and normalize the convolution template signal.

[0007] Step 3: Fast sparse deconvolution restoration. Build a 2D image fast sparse deconvolution restoration model to obtain the deconvolution restored image.

[0008] Image restoration is to use some prior knowledge of the degradation phenomenon to reconstruct or restore the degraded image. The image degradation model is shown in (1):

[0009] G=AF+N (1)

[0010] Where G represents the observed image, F represents the complex original image, A represents the blur operator matrix derived from the convolution kernel, AF represents the convolution of the convolution template and the original image, and N is the noise signal. The image deblurring problem refers to the process of restoring the original image F from the observed image G. The image deblurring problem is transformed into an optimization problem and is expressed as

[0011]

[0012] λ is the adjustment parameter, ‖·‖2 and ‖·‖1 are the L2 norm and L1 norm. To solve the optimization problem (2), we define

[0013]

[0014] Formula (2) can be further rewritten as

[0015]

[0016] The deconvolution operation of the visible image is to solve the real image signal when the observation signal and template signal are known.

[0017] For a given λ>0, choose an initial vector F 0 ∈R mn For k = 0, 1, 2, ..., ..., by F k F is obtained by iterating as follows k+1 :

[0018] ①

[0019] ②F k+1 =A λ (β k )

[0020] where δ k is the k-th iteration step, is the gradient, which is the calculation of A T (AF k -G) and the value corrected by formula (4).

[0021]

[0022] A λ The definition of A λ ([α1,α2,…,α n ] T )=[t λ (α1)),t λ (α2),...,t λ (α n )] T ,t λ (α i)=(|λ-α i |)sgn(α i ). sgn is the sign function.

[0023] The specific process of step 3 is:

[0024] Step 3.1 Initialization. Construct the initial value F of the deconvolution image 0 .

[0025] Step 3.2 Parameter setting: Set the maximum number of iterations K and the initial value of the adjustment parameter λ, and let t be the counter, initially t = 1.

[0026] Step 3.3: Minimize F2(F). This process uses the gradient projection algorithm and the fast spatial discrete convolution algorithm to minimize F2 in the model. The specific implementation includes the following five steps:

[0027] S1 observed image G minus deconvolution image F k The convolution of (k=0,1,2…,K-1) and the convolution kernel T obtains the difference image R, which realizes G-AF k , AF k A vector-matrix representation of two-dimensional finite discrete convolution for ease of analysis and calculation.

[0028] The present invention realizes the discrete convolution algorithm of spatial two-dimensional image. In this step, this function is used to complete the deconvolution image matrix F k Convolution with convolution template T.

[0029] The steps for fast convolution operation of two-dimensional image spatial domain are as follows:

[0030] ① Find the size of the two-dimensional image F, expressed as m×n, and the size of the convolution template T, expressed as (2p+1)×(2q+1), where p and q can be decimals;

[0031] ② Set the initial value of the convolution result ConvF, ConvF is a zero matrix of size (m+2p)×(n+2q);

[0032] ③ Sequentially select elements F(i, j) (1≤i≤m,1≤j≤n) from the two-dimensional image matrix F and multiply them with the template T; then add the product to the modules corresponding to rows i to i+2p and columns j to j+2q in ConvF to obtain a new ConvF;

[0033] ④ Repeat step ③ until all elements in F have completed step ③.

[0034] ⑤Output the convolution result ConvF of F and T.

[0035] The S2 convolution kernel T is flipped 180° and convolved with the difference image R, and then the center of the result is taken and Fk The negative gradient is obtained by the matrix blocks of the same size, that is, A is completed. T (G-AF k )have to

[0036] The spatial domain fast convolution operation method combined in this step completes A T (G-AF k ). The specific steps are:

[0037] ① Find the size of the difference matrix R, expressed as m×n, and the size of the convolution template T, expressed as (2p+1)×(2q+1), where p and q can be decimals;

[0038] ② Use the above-mentioned fast discrete convolution operation method in the spatial domain to perform the convolution operation of T and R to obtain I;

[0039] ③ Extract the matrix blocks from rows 2p+1 to m and columns 2q+1 to n in I, and complete A T (G-AF k )have to

[0040] S3 will be the gradient The elements in the matrix that are greater than 0 correspond to F k The values less than 0 are set to 0, that is,

[0041]

[0042] S4 finds the optimal step length parameter along the echelon descent direction

[0043]

[0044] Where R is the difference calculated in S1 in step 3.3, ‖·‖ F is the Frobenius norm (·) + Defined as (x) + =max{x,0}.

[0045] S5 calculation of β k . in

[0046] Step 3.4: Minimize F1(F) and get F k+1 , and let F k+1 For the new F k The implementation process uses an iterative threshold reduction algorithm to minimize F1, as follows:

[0047] F k+1 =A λ(β k )

[0048] F k =F k+1

[0049] Among them, A λ ([α1,α2,...,α n ] T )=[t λ (α 1) ),t λ (α2),...,t λ (α n )] T , t λ (α i )=(|λ-α i |)sgn(α i ), sgn is the sign function.

[0050] Step 3.5 determines whether t is less than K. If so, t=t+1, and then proceed to step 3.3; otherwise, proceed to step 3.6.

[0051] In step 3.6, the weighted L1 minimization model is further sparsely adjusted and restored to obtain the final deconvolution restored image.

[0052] S1 sets the convergence condition threshold ξ, the maximum number of weighted adjustments p, and the maximum number of sparse inverse convolution deblurring under each weighted adjustment number q, and sets the corresponding counting variables m=1 and n=1.

[0053] S2 in the obtained F k Based on this, update the weight λ(i,j) parameters. ε is a small positive decimal.

[0054] S3 iterative solution updates F k The method is the same as steps 3.3 to 3.4.

[0055] S4 determines whether m is less than q. If so, m=m+1, and then proceeds to step S3 of this process; otherwise, proceeds to step S5.

[0056] S5 determines whether the convergence condition is met. Here, we can set F k and F in step S2 k Whether the 2-norm of the difference is less than a given threshold ξ, but not limited to this. If the condition is met, terminate the iteration and go to step S7, otherwise go to step S6.

[0057] S6 determines whether n is less than p. If it is less than p, n=n+1, and then proceeds to step S2; otherwise, proceeds to step S7.

[0058] S7 sets F = F k , output the restored image F. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1 Flowchart of the fast sparse deconvolution two-dimensional image restoration method provided by the present invention;

[0060] Figure 2 It is the minimized F2(F) module provided by the present invention;

[0061] Figure 3 is the standard test grayscale image Lena in the example of the present invention;

[0062] Figure 4 is the observed image after blurring and adding noise in the example of the present invention;

[0063] Figure 5 is the peak signal-to-noise ratio between the restored image and the true image in the first stage of the embodiment of the present invention;

[0064] Figure 6 is the peak signal-to-noise ratio between the restored image and the true image in the second stage of the embodiment of the present invention;

[0065] Figure 7 It is the sparse deconvolution restored image in the example of the present invention. DETAILED DESCRIPTION

[0066] In order to facilitate those skilled in the art to understand and implement the present invention, the present invention is further described in detail below in conjunction with the accompanying drawings and implementation examples. It should be noted that the implementation examples described herein are only used to illustrate and explain the present invention and are not intended to limit the present invention. Figure 1 This is the process of the fast sparse deconvolution two-dimensional image restoration method provided by the present invention. The specific implementation method is as follows:

[0067] Step 1: Data input. Input the observed image and convolution template. The present invention has no special requirements for the convolution kernel and noise distribution. The following is just an example. To facilitate effect comparison and performance evaluation, we use the standard grayscale test image LENA with a size of 256×256 pixels. Figure 3 As shown, the test image matrix is convolved with the template T and then the noise is added to the observed image G as shown in Figure 4 The template used in this implementation is

[0068]

[0069] Add Gaussian noise with mean 0 and standard deviation 0.02.

[0070] The second step is to linearly transform the value range of the observed image pixels and normalize the convolution template signal.

[0071] Step 3: Fast sparse inverse convolution restoration processing.

[0072] Step 3.1 Initialization. Construct the initial value F of the deconvolution image 0 Construct the initial deconvolution image based on the template size and the observed image. Initial deconvolution image F 0 , whose size is the row and column values of the observed image minus the row and column values of the template plus 1.

[0073] Step 3.2: Parameter Setting. Set the maximum number of iterations K and the initial value of the adjustment parameter λ. Let t be a counter, initially t = 1. In the first phase, the initial value of λ is 0.05 and K = 30.

[0074] Step 3.3: Minimize F2(F). This implementation uses the gradient projection algorithm combined with the fast convolution algorithm disclosed in the present invention to minimize the L2 norm in the model, such as Figure 2 As shown, the algorithm implementation includes the following five steps S1 to S5:

[0075] S1 observed image G minus deconvolution image F k The convolution of (k=0,1,2…,K-1) and the convolution kernel T obtains the difference image R, which realizes G-AF k , AF k A vector-matrix representation of two-dimensional finite discrete convolution for ease of analysis and calculation.

[0076] The present invention realizes the discrete convolution algorithm of spatial two-dimensional image. In this step, this function is used to complete the deconvolution image matrix F k Convolution with convolution template T.

[0077] The steps for fast convolution operation of two-dimensional image spatial domain are as follows:

[0078] ① Find the size of the two-dimensional image F, expressed as m×n, and the size of the convolution template T, expressed as (2p+1)×(2q+1), where p and q can be decimals;

[0079] ② Set the initial value of the convolution result ConvF, ConvF is a zero matrix of size (m+2p)×(n+2q);

[0080] ③ Sequentially select elements F(i, j) (1≤i≤m,1≤j≤n) from the two-dimensional image matrix F and multiply them with the template T; then add the product to the modules corresponding to rows i to i+2p and columns j to j+2q in ConvF to obtain a new ConvF;

[0081] ④ Repeat step ③ until all elements in F have completed step ③.

[0082] ⑤Output the convolution result ConvF of F and T.

[0083] The S2 convolution kernel T is flipped 180° and convolved with the difference image R, and then the center of the result is taken and F k The matrix blocks of the same size are used to find the negative gradient and complete A. T (G-AF k )have to

[0084] The spatial domain fast convolution operation method combined in this step completes A T (G-AF k ). The specific steps are:

[0085] ① Find the size of the difference matrix R, expressed as m×n, and the size of the convolution template T, expressed as (2p+1)×(2q+1), where p and q can be decimals;

[0086] ② Using the spatial domain fast convolution operation method of the present invention to perform convolution operation of T and R to obtain I;

[0087] ③ Extract the matrix blocks from rows 2p+1 to m and columns 2q+1 to n in I, and complete A T (G-AF k )have to

[0088] S3 will be the gradient The elements in the matrix that are greater than 0 correspond to F k The values less than 0 are set to 0, that is,

[0089]

[0090] S4 finds the optimal step length parameter along the echelon descent direction

[0091]

[0092] in,‖·‖ F is the Frobenius norm (·) + Defined as (x) + =max{x,0}.

[0093] S5 calculation of β k . in

[0094] Step 3.4: Minimize F1(F) and get F k+1 , and let F k+1 For the new F kThe implementation process uses an iterative threshold reduction algorithm to minimize F1, as follows:

[0095] F k+1 =A λ (β k )

[0096] F k =F k+1

[0097] Among them, A λ ([α1,α2,...,α n ] T )=[t λ (α 1) ),t λ (α2),...,t λ (α n )] T , t λ (α i )=(|λ-α i |)sgn(α i ), sgn is the sign function.

[0098] Step 3.5 determines whether t is less than K. If it is less than K, then t = t + 1, and then go to step 3.3; otherwise, go to step 3.6. This stage can quickly remove noise and improve the deblurring effect. The peak signal-to-noise ratio of the restored image and the real image in the first stage is Figure 5 shown.

[0099] In step 3.6, the weighted L1 minimization model is further sparsely adjusted and restored to obtain the final result.

[0100] S1 sets the convergence condition threshold ξ, the number of weighted adjustments p, and the number of sparse inverse convolution deblurring operations q for each weighted adjustment, and sets the corresponding counting variables m = 1 and n = 1. In this implementation, p is 10 and q is 30.

[0101] S2 in the obtained F k Based on this, update the weight λ(i,j) parameters. ε is a small positive decimal.

[0102] S3 iterative solution updates F k The method is the same as steps 3.3 to 3.4.

[0103] S4 determines whether m is less than q. If so, m=m+1, and then proceeds to step S3 of this process; otherwise, proceeds to step S5.

[0104] S5 determines whether the convergence condition is met. The convergence condition of this example is F k and F in step S2k Is the 2-norm of the difference less than the given convergence threshold ξ? If the condition is met, terminate the iteration and go to step S7, otherwise go to step S6.

[0105] S6 determines whether n is less than p. If it is less than p, n=n+1, and then proceeds to step S2; otherwise, proceeds to step S7.

[0106] S7 sets F = F k , output the deconvolution restoration image F. Compared with the first stage, this stage improves the L1 minimization model and performs weighted threshold reduction, which can greatly improve the restoration effect. The peak signal-to-noise ratio of the restored image in the second stage to the original image is as follows: Figure 6 As shown. In this example, the deconvolution restored image is as follows Figure 7 shown.

[0107] The above description is only a preferred embodiment of the present invention, but the present invention should not be limited to the contents disclosed in the embodiment and the accompanying drawings. Therefore, any equivalent or modified implementation that does not depart from the spirit disclosed in the present invention shall fall within the scope of protection of the present invention.

Claims

1. A fast sparse deconvolution two-dimensional image restoration method, characterized by: The steps include: Step 1: Data input, read the observed image and convolution template; Step 2: Linearly transform the value range of the observed image pixels and normalize the convolution template signal; The third step is fast sparse inverse convolution restoration processing; construct a two-dimensional image fast sparse inverse convolution method restoration model to obtain the inverse convolution restoration image; let G be the observed image, the original image be F, A represent the blur operator matrix, AF is the convolution of the convolution function and the original image; the image restoration model is expressed as solving the optimization problem: λ is the adjustment parameter, ‖·‖2 and ‖·‖1 are the L2 norm and L1 norm; Define F1(F) = λ‖F‖1 and The image restoration model is rewritten as s.t.F i,j ≥0 The specific process of step 3 is: Step 3.1 Initialization: Construct the initial value F of the deconvolution image 0 ; Step 3.2: Parameter setting: set the maximum number of iterations K and the initial value of the adjustment parameter λ, and let t be a counter with the initial t being 1; Step 3.3 minimize F2(F); Step 3.4: Minimize F1(F) and get F k+1 , and let F k+1 For the new F k , and its implementation process is as follows: F k+1 =A λ (b k ) F k =F k+1 Among them, A λ ([a1,a2,...,a n ] T )=[t λ (a 1) ),t λ (a2),...,t λ (a n )] T ,t λ (a i )=(|λ-a i |)sgn(a i ) Step 3.5: Determine whether t is less than K. If so, t = t + 1, and go to step 3.3; otherwise, go to step 3.

6. In step 3.6, the weighted L1 minimization model is further sparsely adjusted and restored to obtain the final result.

2. The fast sparse deconvolution two-dimensional image restoration method according to claim 1, characterized in that: In step 3.3, the algorithm for minimizing the F2(F) module includes the following five steps: S1 observed image G minus deconvolution image F k Convolution of (k=0,1,2…,K-1) with convolution kernel T to achieve G-AF k Get the difference image R; The S2 convolution kernel T is flipped 180° and convolved with the difference image R, and then the center of the result is taken and F k The matrix blocks of the same size are used to find the negative gradient and complete A. t (G-AF k )have to S3 will be the gradient The elements in the matrix that are greater than 0 correspond to F k The values less than 0 are set to 0, that is, S4 finds the optimal step length parameter along the echelon descent direction in,‖·‖ F is the Frobenius norm (·) + Defined as (x) + =max{x,0} S5 calculation of β k , in 3. The fast sparse deconvolution two-dimensional image restoration method according to claim 1, characterized in that: In step 3.6, S1 sets the convergence condition threshold ξ, the maximum number of weighted adjustments p, and the number of sparse inverse convolution deblurring operations q for each weighted adjustment, and sets the corresponding counting variables m=1 and n=1; S2 in the obtained F k Based on this, update the weight λ(i,j) parameter, ε is a small positive decimal; S3 iterative solution updates F k , the method is the same as steps 3.3 to 3.4; S4 determines whether m is less than q. If so, m=m+1, and then proceeds to step S3; otherwise, proceeds to step S5. S5 determines whether the convergence condition is met. Here, we can set F k and F in step S2 k Whether the 2-norm of the difference is less than a given threshold ξ, but not limited to this; if the condition is met, terminate the iteration and go to step S7, otherwise go to step S6; S6 determines whether n is less than p. If so, n=n+1, and then the process goes to step S2; otherwise, it goes to step S7. S7 sets F = F k , output the restored image F.

4. The fast sparse deconvolution two-dimensional image restoration method according to claim 1, characterized in that: The image restoration process is processed from coarse to fine; the processing process is divided into two major stages: the first stage is the coarse processing stage, which uses a unified adjustment parameter λ to minimize the L1 norm. While restoring the image, it is more important to denoise, which is characterized by steps 3.3 to 3.5; the second stage is the fine processing stage, which uses a weighted L1 minimization model to minimize the L1 norm. Larger weights are used to prevent near-zero terms in the reconstructed signal, while smaller weights can be used to encourage non-zero terms, reducing the blur caused by the unified adjustment parameter λ minimizing the L1 norm and retaining more detail information, which is characterized by step 3.

6.

5. The fast sparse deconvolution two-dimensional image restoration method according to claim 2, characterized in that: Deconvoluted image F in step S1 of step 3.3 k The convolution operation of (k=0,1,2…,K) and the convolution kernel T is performed. This convolution operation adopts the fast convolution operation of the two-dimensional digital image spatial domain. The fast convolution operation of the two-dimensional digital image spatial domain is characterized by: Step 1: Find the size of the two-dimensional image F, expressed as m×n, and the size of the convolution template T, expressed as (2p+1)×(2q+1), where p and q can be decimals; Step 2: Set the initial value of the convolution result ConvF, ConvF is a zero matrix of size (m+2p)×(n+2q); Step 3: Select elements F(i, j) (1≤i≤m,1≤j≤n) from the two-dimensional image matrix F and multiply them with the template T; and add the product to the modules corresponding to rows i to i+2p and columns j to j+2q in ConvF to obtain a new ConvF; Step 4: Repeat step 3 until all elements in F have completed step 3. Step 5 outputs the convolution result ConvF of F and T.

6. The fast sparse deconvolution two-dimensional image restoration method according to claim 2, characterized in that Step S2 of step 3.3 quickly solves the negative gradient; step S2 specifically includes the following process: Step 1: Find the size of the difference matrix R, expressed as m×n, and the size of the convolution template T, expressed as (2p+1)×(2q+1), where p and q can be decimals. Step 2, according to claim 3, performs a spatial convolution operation on T and R to obtain I; Step 3: Extract the matrix blocks with rows 2p+1 to m and columns 2q+1 to n in I, thus completing A. T (G-AF k )have to

Images