Image reconstruction method and system based on primal duality and Anderson acceleration

Abstract

This invention proposes an image reconstruction method and system based on primal-dual and Anderson acceleration. The method comprises the following steps: Step 1: Collect sampled data and establish a signal model; Step 2: Solve for the extreme value after adding the #imgabs0# norm and the TV norm, and the solution provides an estimate of the image x; Step 3: Expand the norm using the primal-dual method and introduce new iterative variables; Step 4: Calculate the estimated value of the variable using the proximal gradient descent method; Step 5: Use the estimated value obtained in Step 4 to perform Anderson acceleration, iterate several times, and obtain the reconstructed image x. Simulation experiments show that the technical solution provided by this invention performs relatively well in terms of estimation accuracy and noise immunity.

Description

Technical Field

[0001] The present invention belongs to the technical field of image information processing, and in particular relates to an image reconstruction method and system based on primal-dual and Anderson acceleration. Background Art

[0002] Image reconstruction has always been a challenging technique in image processing. The goal of image reconstruction is to reconstruct missing information from given information and restore image quality. It is a fundamental and important task in image processing, with widespread applications in various fields such as medicine, military, astronomy, and geophysics. Ensuring high-quality image reconstruction from limited or corrupted (e.g., noisy) measurements (e.g., reducing reconstruction time) has been an active research area.

[0003] In the past few decades, image reconstruction technology has made significant progress, and researchers in this field have proposed many effective methods, such as least squares (LS), maximum likelihood, machine learning, etc. Reconstruction is usually regarded as the solution of a linear inverse problem. Sparse images often appear in nature, for example, in radio astronomy and molecular imaging. Similarly, non-sparse images may have sparse representations in some appropriate domains. For example, an image composed of several fixed-value regions will become sparse after spatial differentiation. As long as the image is sparse in the transform domain or dictionary, compressed sensing (CS) technology can be used to reconstruct the image. There are now several methods to solve the problem of sparse image reconstruction, such as Bayesian learning, least absolute shrinkage and selection operator (LASSO), and fast iterative shrinkage threshold algorithm (FISTA).

[0004] Furthermore, total variation (TV) signal denoising is a commonly used nonlinear filtering method for estimating piecewise constant / smooth signals. TV regularization (penalty) is also widely used in image reconstruction and denoising due to its ability to preserve edges and involves non-quadratic functions of the differences between adjacent pixels. Current image reconstruction methods for solving the inverse problem of total variation regularization include maximization-minimization (MM), primal-dual hybrid gradient (PDHG), and the alternating direction method of multipliers (ADMM).

[0005] In the field of image analysis and processing, there are many typical data that have both sparse features and piecewise constant features. Therefore, when solving the corresponding image reconstruction problem, sparse penalty and TV regularization should be introduced simultaneously to force the solution of the inverse problem with these prior properties. However, the resulting image reconstruction problem involves incorporating multiple non-differentiable regularization terms into the optimization. This increases the complexity of the entire scheme and brings a challenge, namely how to obtain a more efficient method to solve this optimization problem. Although ADMM is a powerful splitting algorithm and is very suitable for this optimization problem by dividing the multi-block separable convex optimization into several subproblems, the subproblems are often difficult to solve effectively due to the lack of closed solutions or the need for large computational costs (such as matrix inversion operations). Summary of the Invention

[0006] In order to solve the problem of reconstructing the original image from a small and noisy received signal in the prior art, the present invention provides an image reconstruction method and system based on primal duality and Anderson acceleration.

[0007] The technical solutions adopted by the present invention to solve the technical problems are as follows:

[0008] An image reconstruction method based on primal-dual and Anderson acceleration comprises the following steps:

[0009] Step 1: Collect sampling data and build a signal model;

[0010] Step 2: Solve the join The extreme values ​​after the norm and TV norm, the solution provides an estimate of the image x;

[0011] Step 3: Use the primal-dual method to expand the norm and introduce new iterative variables;

[0012] Step 4: Use proximal gradient descent to calculate the estimated value of the variable;

[0013] Step 5: Use the estimated value obtained in step 4 to perform Anderson acceleration and iterate several times to obtain the reconstructed image x.

[0014] Preferably, in step 1, least squares (LS) is introduced, and image reconstruction begins with the form of least squares (LS) solution to the problem:

[0015]

[0016] Among them, ||·|| is the Euclidean norm; is the discretization of the integral operator, Measurements corrupted by errors are collected; a signal model is built based on the data y:

[0017] y=Hx0+v

[0018] in, Represents an unknown two-dimensional image with N pixels, whose entries have been reordered in lexicographic order. Collect all the interference; let the image x0 be a sparse piecewise constant, and each item of v is a realization of a uniformly distributed white noise random variable. As a preference, in step 2, solve the problem of adding The extreme value after the norm and TV norm, due to the noise in the observation data y and the serious pathology of H, the solution in image reconstruction is the difference approximation of x0, where represents the Moore-Penrose pseudo-inverse of H; compressed sensing total variation regularization is used, that is:

[0019]

[0020] Where, It is used to enhance signal sparsity norm, ||x|| TV is the TV norm used to maintain edge regularization and is defined as:

[0021]

[0022] in,

[0023] Among them, σ 2 is the variance, Δ ij is a function of the physical or feature space distance between the i-th point and the j-th point. Preferably, in step 3, the primal-dual method is used to expand the norm and introduce a new iteration variable. For any given norm ||·||, its dual norm is defined as ||x|| * =max ||z||≤1 <x,z> ; In addition, there is the following relationship between the original module and the dual module: 1 and Therefore, we can get:

[0024]

[0025]

[0026] Where S={Z=[z1,…,z N ] T |||z i ||2≤1,i=1,…,N} and Z=[z1,…,z N ] T , these two sets S and Q are convex; written in concave-convex saddle point optimization form:

[0027]

[0028] Where x is the two-dimensional image of the original variable and the position to be solved, and Z and h are the dual variables.

[0029] As a preference, in step 4, the estimated values ​​of the variables x, Z, and h are calculated using the proximal gradient descent method; first, a divergence operator div is introduced, therefore have to:

[0030]

[0031] And, given any convex function h, its nearest neighbor operator is defined as:

[0032]

[0033] x (K) ,h (K) ,Z (K) represents an indicator function, expressed as:

[0034]

[0035] Using the proximal gradient descent method, the problem-solving framework is as follows:

[0036]

[0037] Preferably, in step 5, the results of the first m iterations are used as the conditions for the next iteration, and a fixed point iteration is first set:

[0038] (x (k+1) ,h (k+1) ,Z (k+1) )=g(x (k) ,h (k) ,Z (k) )

[0039] Where g(·)={g1(·),g2(·),g3(·)}, then:

[0040]

[0041] let

[0042]

[0043] In the kth iteration as the residual vector, the mean is defined as Among them, z i represents the i-th column vector of Z; the previous m iterations are used to construct a new iteration that converges to the point (x(k+1) ,h (k+1) ,Z (k +1) ), involves finding a coefficient vector Make To solve the constrained optimization problem, we can get Among them, Q k is an m k ×m k A matrix where the (i,j)th element Then, define the Lagrangian:

[0044]

[0045] Taking the derivative of α in L(α,λ) and setting the result to zero, we get Because 1 T α=1, so we get the solution to the constrained optimization problem:

[0046]

[0047] The next iteration of each parameter is obtained as:

[0048]

[0049] After several iterations (preferably 500), the reconstructed image x is obtained.

[0050] The present invention also discloses an image reconstruction system based on primal duality and Anderson acceleration, which includes the following modules based on the above method:

[0051] Signal model building module: collects sampling data and builds signal models;

[0052] Image estimation module: Solve the joining The extreme values ​​after the norm and TV norm, the solution provides an estimate of the image x;

[0053] Norm expansion module: uses the primal-dual method to expand the norm and introduce new iterative variables;

[0054] Variable calculation module: uses the proximal gradient descent method to calculate the estimated value of the variable;

[0055] Image reconstruction module: Use the estimated value to perform Anderson acceleration and iterate several times to obtain the reconstructed image x.

[0056] Compared with the prior art, the technical effects of the present invention are:

[0057] This paper considers compression denoising during signal storage and transmission and proposes a new algorithm for image reconstruction. This method uses a primal-dual approach and Anderson acceleration for image reconstruction, and the resulting technical solution demonstrates excellent estimation accuracy and convergence speed. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] Figure 1 A preferred embodiment of the present invention relates to a signal model for image reconstruction;

[0059] Figure 2 This is a specific flow chart of an image reconstruction method based on primal duality and Anderson acceleration according to a preferred embodiment of the present invention;

[0060] Figure 3 This is a comparison of the reconstruction results of different algorithms on the Set14 dataset;

[0061] Figure 4 This is the relationship between the relative error and the number of iterations when the sampling rate is 100%;

[0062] Figure 5 is the relationship between relative error and sampling rate;

[0063] Figure 6 It is a graph of the relationship between the total iteration time and the sampling rate;

[0064] Figure 7 This is the relationship between relative error and signal-to-noise ratio when the sampling is 100%, 90%, and 80% respectively;

[0065] Figure 8 This is the relationship between relative error and number of iterations when the signal-to-noise ratio is 15dB and 5dB respectively;

[0066] Figure 9 This is the imaging result diagram after 50 iterations without and with the Anderson acceleration method;

[0067] Figure 10 This is a block diagram of an image reconstruction system based on primal-dual and Anderson acceleration according to a preferred embodiment of the present invention. DETAILED DESCRIPTION

[0068] The technical solution of the present invention will be further described in detail below through specific embodiments in conjunction with the accompanying drawings.

[0069] like Figure 1-2 As shown, this embodiment discloses an image reconstruction method based on primal duality and Anderson acceleration, which includes the following steps:

[0070] Step 1: Collect sampling data and build a signal model; the details are as follows:

[0071] The problem solved by least squares (LS) is of the form:

[0072]

[0073] Among them, ||·|| is the Euclidean norm; is the discretization of the integral operator, Measurements corrupted by errors are collected; a signal model is built based on the data y:

[0074] y=Hx0+v

[0075] in, represents an unknown two-dimensional image with N pixels whose entries have been reordered lexicographically. All the interferences are collected; let the image x0 be a sparse piecewise constant and each entry of v be a realization of a uniformly distributed white noise random variable.

[0076] Step 2: Solve the join The extreme values ​​after norm and TV norm are as follows:

[0077] Solve Join The extreme value after the norm and TV norm, due to the noise in the observation data y and the serious pathology of H, the solution in image reconstruction is the difference approximation of x0, where represents the Moore-Penrose pseudo-inverse of H; compressed sensing total variation regularization is used, that is:

[0078]

[0079] Where, It is used to enhance signal sparsity norm, ||x|| TV is the TV norm used to maintain edge regularization and is defined as:

[0080]

[0081] in,

[0082] Among them, σ 2 is the variance, Δ ij is a function of the physical or feature space distance between the i-th point and the j-th point.

[0083] Step 3: Use the primal-dual method to expand the norm and introduce new iteration variables, as follows:

[0084] The primal-dual method is used to expand the norm and introduce new iteration variables. For any given norm ||·||, its dual norm is defined as ||x|| * =max ||z||≤1 <x,z> ; In addition, there is the following relationship between the original module and the dual module: and Therefore, we get:

[0085]

[0086] Where S={Z=[z1,…,z N ] T |||zi||2≤1,i=1,…,N} and Z=[z1,…,z N ] T , these two sets S and Q are convex; written in concave-convex saddle point optimization form:

[0087]

[0088] Where x is the original variable and the unknown two-dimensional image to be solved, and Z and h are the dual variables.

[0089] Step 4: Use the proximal gradient descent method to calculate the estimated value of each variable, as follows:

[0090] The estimated values ​​of x, Z, and h variables are calculated using the proximal gradient descent method; a divergence operator div is first introduced. therefore We can get:

[0091]

[0092] And, given any convex function h, its nearest neighbor operator is defined as:

[0093]

[0094] x (K) ,h (K) ,Z (K) represents an indicator function, expressed as:

[0095]

[0096] Using the proximal gradient descent method, the problem-solving framework is as follows:

[0097]

[0098] Step 5: Use the estimated value obtained in step 4 to perform Anderson acceleration, continue multiple iterations, update the parameter x, and obtain the reconstructed image x; the details are as follows:

[0099] Use the results of the first m iterations as the conditions for the next iteration and set a fixed point iteration first:

[0100] (x (k+1) ,h (k+1) ,Z (k+1) )=g(x (k) ,h (k) ,Z (k) )

[0101] Where g(·)={g1(·),g2(·),g3(·)}, then:

[0102]

[0103] let

[0104]

[0105] In the kth iteration as the residual vector, the mean is defined as Among them, z i represents the i-th column vector of Z; the previous m iterations are used to construct a new iteration that converges to the point (x (k+1) ,h (k+1) ,Z (k +1) ), involves finding a coefficient vector Make To solve the constrained optimization problem, we have Among them, Q k is an m k ×m k A matrix where the (i,j)th element Then, define the Lagrangian:

[0106]

[0107] Taking the derivative of α in L(α,λ) and setting the result to zero, we get Because 1 T α=1, so the solution to the constrained optimization problem is obtained:

[0108]

[0109] The next iteration of each parameter is obtained as:

[0110]

[0111] After 500 iterations, the reconstructed image x is obtained.

[0112] This paper evaluates the performance of the proposed image reconstruction based on the primal-dual method (CTPD) by reconstructing images and numerical results, and compares it with other existing methods, such as the LS method, -Regularization (LASSO), TV-regularization, ADMM-based -TV composite regularization method and FISTA-based -Regularization algorithm.

[0113] Figure 3 The reconstruction performance of these algorithms is shown in the case of a signal-to-noise ratio of 15dB in 500 iterations. From left to right, they are Original, LS, LASSO, TV, ADMM, FISTA, and CTPD. The signal-to-noise ratio is defined as SNR: = 10log(P s / P n In this example, λ1 = 0.02, λ2 = 0.02, ν = 0.05, and μ = 0.1 are set, with H being the sampling matrix. Image reconstruction results of various algorithms are shown when the sampling rate is 90% (M / N = 90%). The imaging quality of the CTPD algorithm is superior to that of LS, LASSO, and FISTA, and the visual effect is similar to that of TV and ADMM. Therefore, it can be concluded that the CTPD algorithm has good image reconstruction performance under incomplete data sampling.

[0114] Figure 4 Shows the relative error at different iteration numbers Among them, x and Represent the target model and the reconstructed image respectively. It can be seen that after convergence, CTPD has a lower relative error.

[0115] Figure 5 and Figure 6 The error rates and imaging times of the imaging results obtained by LS, LASSO, TV, ADMM, FISTA and CTPD under different data ratios (sampling rates) are shown respectively. Figure 5 It can be seen that when the data ratio is 40% or 50%, the relative error of ADMM is about 1dB lower than that of CTPD. In other cases, CTPD has a lower error rate. Since the matrix dimension of w is large when CTPD calculates the TV norm, the calculation time is longer, such as Figure 6 shown.

[0116] Figure 7 The relative error of the algorithm at different signal-to-noise ratios (SNRs) is shown for sampling rates of 100%, 90%, and 80%. It can be seen that when the SNR is greater than 5dB, the error rate of the CTPD algorithm is generally low. Furthermore, the advantages of the method of the present invention become increasingly significant as the sampling rate decreases.

[0117] In order to restore the original signal with fewer iterations, the present invention adds the Anderson accelerated method (A-CTPD) on the basis of solving the problem using the primal-dual method. Figure 8 The relative errors of different iteration times are given when the signal-to-noise ratio is 15dB and 5dB respectively. The experimental results show that the A-CTPD algorithm not only has a fast convergence speed, but also has the same image reconstruction accuracy as the CTPD algorithm.

[0118] Figure 9 The following are the imaging results of the CTPD and A-CTPD algorithms after 50 iterations, with parameters set to ν = 0.01 and μ = 0.1. It can be seen that with the same number of iterations, A-CTPD exhibits better imaging quality than CTPD.

[0119] In summary, the technical solution proposed in the present invention adopts a sparse and total variation composite regularization method to minimize multiple cost functions. This minimization problem does not directly find its optimal solution, but first converts it into a convex-concave saddle point optimization problem, and solves it through a primal-dual iterative algorithm with several simple sub-steps, and the solution complexity of each sub-step is low. In addition, in order to speed up the convergence speed, an acceleration algorithm is proposed on this basis. In this algorithm, the primal-dual iteration is restated as a fixed-point iteration, so the Anderson acceleration technology can be used to converge faster without sacrificing the image reconstruction accuracy. Through simulation experiments, it can be obtained that the technical solution provided by the present invention has relatively excellent performance in both estimation accuracy and noise resistance.

[0120] like Figure 10 As shown, this embodiment discloses an image reconstruction system based on primal duality and Anderson acceleration, based on the above method embodiment, including the following modules:

[0121] Signal model building module: collects sampling data and builds signal models;

[0122] Image estimation module: Solve the joining The extreme values ​​after the norm and TV norm, the solution provides an estimate of the image x;

[0123] Norm expansion module: uses the primal-dual method to expand the norm and introduce new iterative variables;

[0124] Variable calculation module: uses the proximal gradient descent method to calculate the estimated value of the variable;

[0125] Image reconstruction module: Use the estimated value to perform Anderson acceleration and iterate several times to obtain the reconstructed image x.

[0126] For other contents of this embodiment, please refer to the above method embodiment.

[0127] It should be understood that the above description of the preferred embodiment is relatively detailed and cannot be regarded as limiting the scope of protection of the patent of the present invention. Under the guidance of the present invention, ordinary technicians in this field can make substitutions or modifications without departing from the scope of protection of the claims of the present invention, which fall within the scope of protection of the present invention. The scope of protection requested by the present invention shall be based on the attached claims.

Claims

1. An image reconstruction method based on primal-dual and Anderson acceleration, characterized in that: The steps include: Step 1: Collect sampling data and build a signal model; Step 2: Find the extreme value after adding the l1 norm and TV norm. The solution provides an estimate of the image x. Step 3: Use the primal-dual method to expand the norm and introduce new iterative variables; Step 4: Use proximal gradient descent to calculate the estimated value of the variable; Step 5: Use the estimated value obtained in step 4 to perform Anderson acceleration and iterate several times to obtain the reconstructed image x; In step 1, the least squares LS problem is solved in the form of: Among them, ||·|| is the Euclidean norm; is the discretization of the integral operator, Measurements corrupted by errors are collected; a signal model is built based on the data y: y=Hx0+v in, Represents an unknown two-dimensional image with N pixels, the entries of which have been reordered in lexicographic order. All interferences are collected; let the image x0 be a sparse piecewise constant, and each item of v is a realization of a uniformly distributed white noise random variable; In step 2, the solution in image reconstruction is the difference approximation of x0, where represents the Moore-Penrose pseudo-inverse of H; compressed sensing total variation regularization is used, that is: Where, is the l1 norm used to enhance signal sparsity, ||x|| TV is the TV norm used to maintain edge regularization and is defined as: in, Among them, σ 2 is the variance, Δ ij is a function of the physical or feature space distance between the i-th point and the j-th point; In step 3, for any given norm ||·||, its dual norm is defined as ||x|| * =max ||z||≤1 <x,z> ; In addition, there is the following relationship between the original module and the dual module: ||x|| 1* =||x|| ∞ ,||x|| ∞* =||x||1 and ||x|| 2* =||x||2; therefore, we get: Where S={Z=[z1,…,z N ] T |||z i ||2≤1,i=1,…,N} and Z=[z1,…,z N ] T , these two sets S and Q are convex; written in concave-convex saddle point optimization form: Where x is the original variable and the unknown two-dimensional image to be solved, and Z and h are the dual variables.

2. The image reconstruction method based on primal-dual and Anderson acceleration according to claim 1, characterized in that: In step 4, the estimated values ​​of x, Z, and h variables are calculated using the proximal gradient descent method; First introduce a divergence operator div, therefore have to: And, given any convex function h, its nearest neighbor operator is defined as: x (K) ,h (K) ,Z (K) represents an indicator function, expressed as: Using the proximal gradient descent method, the problem-solving framework is as follows:

3. The image reconstruction method based on primal-dual and Anderson acceleration according to claim 2, characterized in that: In step 5, the results of the first m iterations are used as the conditions for the next iteration. First, a fixed point iteration is set: (x (k+1) ,h (k+1) ,Z (k+1) )=g(x (k) ,h (k) ,Z (k) ) Where g(·)={g1(·),g2(·),g3(·)}, then: let In the kth iteration as the residual vector, the mean is defined as Among them, z i represents the i-th column vector of Z; the previous m iterations are used to construct a new iteration that converges to the point (x (k+1) ,h (k+1) ,Z (k+1) ), involves finding a coefficient vector Make To solve the constrained optimization problem, we have Among them, Q k is an m k ×m k A matrix where the (i,j)th element Then, define the Lagrangian: Taking the derivative of α in L(α,λ) and setting the result to zero, we get Because 1 T α=1, so the solution to the constrained optimization problem is obtained: The next iteration of each parameter is obtained as: After iteration, the reconstructed image x is obtained.

4. An image reconstruction system based on primal-dual and Anderson acceleration, based on the method according to any one of claims 1 to 3, characterized in that: Includes the following modules: Signal model building module: collects sampling data and builds signal models; Image estimation module: solves the extreme value after adding the l1 norm and TV norm, and the solution provides an estimate of the image x; Norm expansion module: uses the primal-dual method to expand the norm and introduce new iterative variables; Variable calculation module: uses the proximal gradient descent method to calculate the estimated value of the variable; Image reconstruction module: Use the estimated value to perform Anderson acceleration and iterate several times to obtain the reconstructed image x.

Images