Fast robust sparse Bayesian image reconstruction method based on generalized approximate message passing

Abstract

The invention discloses a fast robust sparse Bayesian image reconstruction method based on generalized approximate message passing. The method comprises the following steps: S1, selecting an image data set; s2, setting the pixel sizes of all the images; s3, performing sparse preprocessing on each image; s4, establishing a sparse Bayesian model; s5, calculating a likelihood function, prior distribution of a sparse coefficient x and a marginal likelihood function; s6, calculating the posterior density of x according to the Bayesian theorem; s7, calculating a cost function; s8, iteratively updating the hyper-parameters by adopting a block coordinate descent method; s9, introducing a generalized approximate message passing algorithm to update a mean value and a variance; s10, introducing a damping factor to enhance convergence; s11, judging whether an iteration termination condition is met or not, if yes, outputting posteriori estimation of x, and if not, executing the step S8; and S12, performing inverse wavelet transform on the output x to obtain a reconstructed image. According to the method, reconstruction can be completed in a short time, the reconstruction error is low, and better reconstruction precision and reconstruction efficiency are achieved.

Description

Technical Field

[0001] The invention belongs to the technical field of image reconstruction, and in particular relates to a fast robust sparse Bayesian image reconstruction method based on generalized approximate message passing. Background Art

[0002] In recent years, machine learning algorithms based on sparse representation (SR) have developed rapidly in various fields. The purpose of sparse representation is to represent the main information of the signal with as few non-zero coefficients as possible, so that the characteristic information of the data or signal can be effectively extracted, making it easier for people to obtain important characteristic information of the image, thereby solving the problem of excessive computation caused by high-dimensional feature space and the problem that existing image feature extraction algorithms do not adequately describe image features. Among them, as the inverse problem of SR, sparse signal recovery and sparse reconstruction have received widespread attention in recent years.

[0003] The essence of sparse reconstruction is to estimate and reconstruct sparse data or signals through mathematical models. Common sparse reconstruction methods can be roughly divided into the following three categories: greedy pursuit algorithm, convex relaxation algorithm and Bayesian algorithm. Among them, the greedy pursuit algorithm mainly reconstructs sparse signals by cyclically iteratively selecting the basis vector that contributes the most to the signal. The algorithm complexity is low but the reconstruction performance is poor. The convex relaxation algorithm solves the optimal approximation of the signal through convex optimization and can obtain a unique solution, but the parameters need to be set manually. The Bayesian algorithm is based on the Bayesian theorem and completes approximate inference and solves problems through probability distribution. As a typical algorithm theory, the Bayesian algorithm has higher robustness and reconstruction accuracy than the other two algorithms.

[0004] Image reconstruction is one of the application fields of sparse reconstruction, and its development in the field of machine vision also plays a vital role. Improving the stability, accuracy and speed of the reconstruction algorithm has always been a research hotspot and difficulty in the field of sparse reconstruction. Among them, the sparse Bayesian Learning algorithm (SBL) is one of the most commonly used methods in image reconstruction. Based on the Bayesian idea, it performs probabilistic modeling on the unknown parameters in the model, and then obtains a sparse solution through Bayesian reasoning, and solves the posterior distribution under the assumed Gaussian prior distribution through the expectation maximization (EM) algorithm.

[0005] The classic SBL algorithm uses prior information such as the sparsity or correlation of the solution to improve the signal recovery performance through the prior model, and is insensitive to the coherence of the measurement matrix, which has advantages in finding the global optimal solution. However, the SBL algorithm has the following problems when reconstructing images: 1) High computational complexity: The EM algorithm needs to calculate the inverse of the high-dimensional matrix at each iteration, and the computational complexity is high on large-scale data sets; 2) Slow convergence: The EM algorithm converges slowly in high-dimensional data space or when the amount of data is large, and it is easy to fall into the local optimal solution.

[0006] The generalized approximate message passing algorithm (GAMP) proposed in recent years completes inference by passing information between nodes of iterative probabilistic graph models, and applies quadratic approximation and Taylor technique approximation to cyclic belief propagation, which can reduce computational complexity and improve reconstruction efficiency. Based on this, the present invention proposes a fast robust sparse Bayesian image reconstruction method based on generalized approximate message passing, called the GAMP-FRSBL method, which can effectively improve the image reconstruction accuracy. Summary of the invention

[0007] In view of the above-mentioned status of the prior art, the present invention proposes a fast and robust sparse Bayesian image reconstruction method based on generalized approximate message passing. The present invention first introduces a generalized approximate message passing algorithm on the basis of SBL to replace the E step of the classic EM algorithm, thereby avoiding the matrix inversion problem and greatly reducing the computational complexity of the model; and introduces a damping factor to enhance the convergence of the model. Finally, combined with the convex optimization strategy, the hyperparameters are iteratively updated by the block coordinate descent (BCD) method to improve the reconstruction accuracy of the model. The present invention improves the image reconstruction accuracy and convergence speed by using the generalized approximate message passing algorithm and the block coordinate descent method. Compared with the existing image reconstruction technology, the image reconstruction using the technical solution of the present invention has significant improvements in reconstruction accuracy and efficiency, and can well meet the application requirements in actual scenarios.

[0008] The present invention adopts the following technical scheme:

[0009] A fast and robust sparse Bayesian image reconstruction method based on generalized approximate message passing comprises the following steps:

[0010] S1, select image dataset;

[0011] S2, set the pixel size of all images;

[0012] S3, performing sparse preprocessing on each image;

[0013] S4, establish a sparse Bayesian model;

[0014] S5, calculating the likelihood function, the prior distribution of the sparse coefficient x, and the marginal likelihood function;

[0015] S6, calculating the posterior density of the sparse coefficient x by Bayes' theorem according to the likelihood function, prior distribution and marginal likelihood function calculated in step S5;

[0016] S7, calculating the cost function;

[0017] S8, using block coordinate descent (BCD) to iteratively update the hyperparameters;

[0018] S9, introduce the generalized approximate message passing (GAMP) algorithm to update the mean and variance;

[0019] S10, introducing damping factor to enhance convergence;

[0020] S11, judging whether the iteration termination condition is met, if so, outputting the estimation of the sparse coefficient x; if not, continuing to execute step S8 until the iteration termination condition is met;

[0021] S12. Perform inverse wavelet transform on the output sparse coefficient x to obtain a reconstructed image.

[0022] The reconstruction performance can be evaluated by five indicators: average running time (TIME), reconstruction error (ERROR), structural similarity (SSIM), peak signal-to-noise ratio (PSNR) and average number of iterations (Mean Iterations).

[0023] Preferably, in step S2, the pixel size of all images is set to 256×256 and processed into .mat format files to store the feature data of the images.

[0024] Preferably, in step S3, a "symmlet,8" wavelet transform with a coarsest scale of 3 and a finest scale of 7 is added to each image for sparse preprocessing, and a multi-scale compressed sensing (CS) scheme is used for image reconstruction. The wavelet coefficients to be estimated are set to: n = 16320, and the number of compressed samples is set to: m = 10608.

[0025] Preferably, step S4, establish a sparse Bayesian (SBL) mathematical model:

[0026] y=Φx+e

[0027] Where y∈R m is the observed data, x∈R n is the sparse coefficient, Φ∈R m×nis a known dictionary matrix, e∈R m is the noise vector. According to Bayes’ theorem, the estimate of x is obtained by maximizing the posterior distribution p(x|y; α, λ):

[0028]

[0029] Preferably, in step S5, the likelihood function p(y|x; λ) is calculated according to the posterior distribution expression in step S4:

[0030] Assume that the noise e has a mean of 0 and a variance of λ -1 Gaussian distribution of I. Where I represents the identity matrix, that is:

[0031] The hyperparameter λ follows the Gamma prior: Among them, the shape parameter c>0 and the rate parameter d>0. is the Gamma function.

[0032] Therefore, the Gaussian likelihood function expression of the observed data y can be obtained as:

[0033]

[0034] Preferably, in step S5, the prior distribution p(x; α) is calculated according to the posterior distribution expression in step S4:

[0035]

[0036] in, is a diagonal matrix, α i is control x i Hyperparameter of sparsity, i=1,2,…,n. Hyperparameter α i Obey the Gamma prior distribution:

[0037]

[0038] Among them, the shape parameter a>0 and the rate parameter b>0. is the Gamma function.

[0039] Preferably, step S5 calculates the marginal likelihood function p(y; α, λ) according to the posterior distribution expression in step S4:

[0040]

[0041] According to the properties of Gaussian function It can be obtained that the marginal likelihood function p(y; α, λ) approximately obeys the mean of 0, and the covariance matrix is Gaussian distribution of . Among them, the matrix F is a matrix, g is a vector, and C is a constant.

[0042] Preferably, in step S6, according to Bayes’ theorem, the convolution result of two Gaussian functions is still a Gaussian function. Therefore, it can be calculated that the posterior density p(x|y; α, λ) of x also obeys the Gaussian distribution, and its posterior mean is and the posterior covariance It can be expressed as:

[0043]

[0044] Preferably, in step S7, the cost function is calculated according to the posterior density obtained in step S6. Specifically, the cost function is as follows:

[0045] λ and Λ in step S6 α These are all hyperparameters to be estimated. Type II maximum likelihood estimation can be used to update α and λ by maximizing the joint probability density function p(y, α, λ):

[0046]

[0047] Maximizing the joint probability density function is equivalent to minimizing -2ln(α,λ). Let σ 2 =λ -1 , Λ γ =diag(γ i )=Λ α ; Ignore the constant term and get the cost function L that combines x and hyperparameters x (x,γ,σ 2 )for:

[0048]

[0049] Preferably, step S8, by iteratively minimizing the cost function L x (x,γ,σ 2 ), and use block coordinate descent (BCD) to adjust the hyperparameters γ and σ 2 To update:

[0050]

[0051] Preferably, in step S9, a generalized approximate message passing (GAMP) algorithm is introduced to update the mean in a low-complexity manner. and variance

[0052] For the prior p(x; α) of x and the likelihood function p(y|x; λ), GAMP uses Taylor expansion and quadratic approximation to approximate the maximum a posteriori estimate (MAP) of x:

[0053]

[0054] Among them, q represents the system input vector (prior information), and y represents the system output vector (posterior information). Therefore, the input function and output function are:

[0055]

[0056] Among them, the variable and Respectively Approximate estimates and noise-free models of The approximate values ​​of and Let S = |Φ| 2 , S represents the square of the component magnitude of the dictionary matrix Φ. is the mean vector, Represents the covariance matrix The vector consisting of the diagonal elements of . and represents the mean and variance of the intermediate variable s. Therefore, the update of the input function and the output function can be obtained as:

[0057]

[0058] The iterative update of the mean and variance can be expressed as:

[0059] Preferably, in step S10, a damping factor is introduced to enhance the convergence of the model. The formula for solving the damping coefficient is as follows:

[0060]

[0061] Where A0, B0 and C0 are constant terms. Damping factor θ s and θ x is the root of the above formula, θ s Used to control θ x Used to control

[0062] Preferably, step S11 determines whether x satisfies the iterative convergence condition, and if so, outputs an estimate of x: The iterative convergence condition is expressed as:

[0063]

[0064] Preferably, in step S12, the output x is subjected to an inverse wavelet transform to obtain a reconstructed image. The calculation formula for the reconstruction error is: in, represents the reconstructed image, and P represents the real image. Finally, the reconstruction performance is evaluated by five indicators: average running time (TIME), reconstruction error (ERROR), structural similarity (SSIM), peak signal-to-noise ratio (PSNR), and average number of iterations (MeanIterations).

[0065] The beneficial effects of the present invention are:

[0066] (1) The present invention avoids matrix inversion by introducing the generalized approximate message passing algorithm (GAMP), saves running time, and improves the reconstruction efficiency of the model; and adopts an appropriate damping coefficient to improve the convergence of the model; finally, combined with the convex optimization strategy, the block coordinate descent method (BCD) is used to iteratively update the hyperparameters, thereby improving the robustness of the model and the reconstruction accuracy.

[0067] (2) Compared with the image reconstruction method in the prior art, the present invention can complete the reconstruction in a shorter time with lower reconstruction error, and has better reconstruction accuracy and reconstruction efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0068] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0069] Figure 1 It is a model structure diagram of a fast and robust sparse Bayesian image reconstruction method based on generalized approximate message passing in a preferred embodiment of the present invention;

[0070] Figure 2 It is an algorithm flow chart related to a fast and robust sparse Bayesian image reconstruction method based on generalized approximate message passing in a preferred embodiment of the present invention;

[0071] Figure 3 is a reconstructed image sample of the image data set used in the present invention;

[0072] Figure 4 It is a comparison chart of the reconstruction performance of the image data set used in the present invention. DETAILED DESCRIPTION

[0073] The following describes the implementation of the present invention through specific embodiments, and those skilled in the art can easily understand other advantages and effects of the present invention from the contents disclosed in this specification. The present invention can also be implemented or applied through other different specific implementations, and the details in this specification can also be modified or changed in various ways based on different viewpoints and applications without departing from the spirit of the present invention.

[0074] The preferred embodiment of the present invention discloses a fast robust sparse Bayesian image reconstruction method based on generalized approximate message passing, referring to Figure 1 As shown in the figure, it is the model structure diagram of this method. Figure 2 , which is the algorithm flow chart involved in this method. Figure 3 , which is a reconstructed image sample of the image data set used in the present invention. All images are processed according to the process of the method of this embodiment to illustrate in detail the effect of the method of the present invention on improving the reconstruction accuracy.

[0075] This embodiment provides a fast and robust sparse Bayesian image reconstruction method based on generalized approximate message passing, comprising the following steps:

[0076] S1. Select a public image dataset.

[0077] S2. Set the pixel size of all images to 256×256 and process them into .mat format files to save the feature data of the images.

[0078] S3. Add wavelet transform to each image for sparse preprocessing, and use multi-scale compressed sensing (CS) scheme to reconstruct the image.

[0079] S4. Establish a mathematical model of sparse Bayesian (SBL): y = Φx + e. Where y∈R m is the observed data, x∈R n is the sparse coefficient, Φ∈R m×n is a known dictionary matrix, e∈R m is the noise vector. According to Bayes’ theorem, the estimate of x is obtained by maximizing the posterior distribution p(x|y; α, λ):

[0080]

[0081] S5. Calculate the likelihood function according to the posterior distribution expression in step S4:

[0082]

[0083] Among them, the noise e has a mean of 0 and a variance of λ -1Gaussian distribution of I, I represents the identity matrix. The hyperparameter λ follows the Gamma prior:

[0084]

[0085] Among them, the shape parameter c>0 and the rate parameter d>0. is the Gamma function.

[0086] According to the posterior distribution expression in step S4, calculate the prior distribution:

[0087]

[0088] in, is a diagonal matrix, α i is control x i Hyperparameter of sparsity, i=1,2,…,n. Hyperparameter α i Obey the Gamma prior distribution:

[0089]

[0090] Among them, the shape parameter a>0 and the rate parameter b>0. is the Gamma function.

[0091] According to the posterior distribution expression in step S4, the marginal likelihood function p(y; α, λ) is calculated:

[0092]

[0093] According to the properties of Gaussian function The marginal likelihood function p(y; α, λ) is approximately 0, and the covariance matrix is Gaussian distribution of . Among them, the matrix F is a matrix, g is a vector, and C is a constant.

[0094] S6. According to Bayes' theorem, we can calculate that the posterior density of x p(x|y; α, λ) also obeys Gaussian distribution, and its posterior mean is and the posterior covariance It can be expressed as:

[0095]

[0096]

[0097] S7, calculate the cost function according to the posterior density obtained in S6. The details are as follows:

[0098] λ and Λ in step S6 αThese are all hyperparameters to be estimated. Type II maximum likelihood estimation can be used to update α and λ by maximizing the joint probability density function p(y, α, λ):

[0099]

[0100] Maximizing the joint probability density function is equivalent to minimizing -2ln(α,λ); let σ 2 =λ -1 , Λ γ =diag(γ i )=Λ α . Get the cost function L x (x,γ,σ 2 )for:

[0101]

[0102] S8, iteratively minimize the cost function L x (x,γ,σ 2 ), and use block coordinate descent (BCD) to adjust the hyperparameters γ and σ 2 To update:

[0103]

[0104] S9. The generalized approximate message passing (GAMP) algorithm is introduced to update the mean and variance in a low-complexity manner.

[0105] For the prior p(x; α) of x and the likelihood function p(y|x; λ), GAMP uses Taylor expansion and quadratic approximation to approximate the maximum a posteriori estimate (MAP) of x:

[0106]

[0107] Among them, q represents the system input vector (prior information), and y represents the system output vector (posterior information). Therefore, the input function and output function are:

[0108]

[0109] Among them, the variable and Respectively Approximate estimates and noise-free models of The approximate values ​​of are and Let S = |Φ| 2 , S represents the square of the component magnitude of the dictionary matrix Φ. is the mean vector, Represents the covariance matrix The vector consisting of the diagonal elements of . and represents the mean and variance of the intermediate variable s. Therefore, the update of the input function and the output function can be obtained as:

[0110]

[0111] The iterative update of the mean and variance can be expressed as:

[0112] S10. Introduce the damping factor to enhance the convergence of the model. The formula for solving the damping coefficient is as follows:

[0113]

[0114] Where A0, B0 and C0 are constant terms. Damping factor θ s and θ x is the root of the above formula, θ s For control θ x For control

[0115] S11. Determine whether x satisfies the iterative convergence condition. If so, output the estimate of x: If not, execute step S8; the iterative convergence condition is expressed as:

[0116]

[0117] S12, perform inverse wavelet transform on the output x to obtain a reconstructed image. The calculation formula for the reconstruction error is: in, represents the reconstructed image, and P represents the real image.

[0118] Finally, the reconstruction performance is evaluated by five indicators: average running time (TIME), reconstruction error (ERROR), structural similarity (SSIM), peak signal-to-noise ratio (PSNR) and average number of iterations (Mean Iterations).

[0119] The following is a further detailed description of each step of the above embodiment in conjunction with specific data.

[0120] In step S1, the obtained data sets are specifically: DOTA optical remote sensing image data set, COCO image data set and UC Merced Land Use (UCM) remote sensing image data set. 50 images (a total of 150 images) are randomly selected from each data set for experiment.

[0121] In step S2, the pixel size of all images is set to 256×256, 50 images are randomly selected from each data set and numbered in sequence, and read into .mat files to save the feature data of the images.

[0122] Step S3, add the wavelet transform of "symmlet,8" with the coarsest scale of 3 and the finest scale of 7 to each image for sparse preprocessing, and use the multi-scale compressed sensing (CS) scheme to reconstruct the image. The wavelet coefficients to be estimated are set to: n = 16320, and the number of compressed samples is set to: m = 10608.

[0123] Step S4, based on the sparse representation theory (SR) and Bayesian theorem, establish a sparse Bayesian (SBL) mathematical model: y = Φx + e. Where y∈R m is the observed data, x∈R n is the sparse coefficient, Φ∈R m×n is a known overcomplete dictionary matrix, e∈R m is the noise vector. According to Bayes’ theorem, the estimate of x is obtained by maximizing the posterior distribution p(x|y; α, λ):

[0124]

[0125] Step S5, calculate the likelihood function p(y|x; λ) according to the posterior distribution expression in step S4. A hierarchical multivariate model is used to represent noise under the SBL framework. The noise vector e has a mean of 0 and a variance of λ -1 Gaussian distribution of I. Where I represents the identity matrix, that is:

[0126] Assume that the hyperparameter λ obeys the Gamma prior:

[0127]

[0128] Among them, the shape parameter c>0 and the rate parameter d>0. is the Gamma function.

[0129] Therefore, the Gaussian likelihood function expression of the observed data y can be obtained as:

[0130]

[0131] According to the posterior distribution expression in step S4, calculate the prior distribution p(x; α). Assume that each element in the sparse coefficient x is independent and identically distributed. Add a prior density to x to achieve sparsity. According to the SBL theory, assume that x has a mean of 0 and a variance of Gaussian prior distribution of .

[0132]

[0133] in, is a diagonal matrix, α i is control x i The hyperparameter of sparsity, i = 1, 2, ..., n. Let the hyperparameter α i Obey the Gamma prior distribution:

[0134]

[0135] Among them, the shape parameter a>0 and the rate parameter b>0. is the Gamma function.

[0136] According to the posterior distribution expression in step S4, the marginal likelihood function p(y; α, λ) is calculated. The marginal likelihood function p(y; α, λ) is calculated by the total probability formula:

[0137]

[0138] According to the properties of Gaussian function It can be obtained that the marginal likelihood function p(y; α, λ) approximately obeys the mean of 0, and the covariance matrix is Gaussian distribution of . Among them, the matrix F is a matrix, g is a vector, and C is a constant.

[0139] Step S6, according to Bayes’ theorem, the convolution result of two Gaussian functions is still a Gaussian function, so we can calculate that the posterior density p(x|y; α, λ) of x also obeys the Gaussian distribution, and its posterior mean is and the posterior covariance It can be expressed as:

[0140]

[0141] Step S7, calculate the cost function according to the posterior density obtained in step S6. The details are as follows: λ and Λ in step S6 α are all hyperparameters to be estimated. Type II maximum likelihood estimation can be used to update α and λ by maximizing the joint probability density function p(y, α, α). Therefore, the solution formulas for α and λ are as follows:

[0142]

[0143] Add -2ln(·) transformation to (α,α), so maximizing the joint probability density function p(y,α,α) is equivalent to minimizing -2 ln(α,λ). Let σ 2 =λ -1 , Λγ =diag(γ i )=Λ α Ignoring the constant term, we can get the cost function L that combines x and hyperparameters x (x,γ,σ 2 )for:

[0144]

[0145] Step S8, minimizing the cost function L by iteration x (x,γ,σ 2 ), we can get γ and σ 2 The update is:

[0146]

[0147] It can be seen that the cost function L x (x,γ,σ 2 ) has a partially decomposable structure, so the block coordinate descent (BCD) method is used to update the hyperparameters.

[0148]

[0149] For the cost function L x (x,γ,σ 2 ) Find the minimum value at zero gradient, and we can get γ and σ 2 The update is:

[0150]

[0151] Step S9, introduce the generalized approximate message passing (GAMP) algorithm to update the mean and variance in a low-complexity manner. For the prior p(x; α) of x and the likelihood function p(y|x; λ), GAMP uses Taylor expansion and quadratic approximation to approximate the maximum a posteriori estimate (MAP) of x:

[0152]

[0153] Among them, q represents the system input vector (prior information), and y represents the system output vector (posterior information). Therefore, the input function and output function are:

[0154]

[0155] Among them, the variable and Respectively Approximate estimates and noise-free models of The approximate values ​​of and Let S = |Φ|2 , S represents the square of the component magnitude of the dictionary matrix Φ. is the mean vector, Represents the covariance matrix The vector consisting of the diagonal elements of . and represents the mean and variance of the intermediate variable s. Therefore, the update of the input function and the output function can be obtained as:

[0156]

[0157] The iterative update of the mean and variance can be expressed as:

[0158] Step S10, introducing a damping factor to enhance the convergence of the model. The damping coefficient solution formula is as follows:

[0159]

[0160] Where A0, B0 and C0 are constant terms. Damping factor θ s and θ x is the root of the above formula, θ s For control θ x For control

[0161] Step S11, determine the mean in step S9 Whether the iterative convergence condition is met, if so, output the estimate of x: If not, execute step S8; the iterative convergence condition is expressed as:

[0162]

[0163] Step S12, performing inverse wavelet transform on the output x obtained in step S11 to obtain a reconstructed image. The calculation formula for the reconstruction error is: in, represents the reconstructed image, and P represents the real image.

[0164] Finally, the reconstruction performance of the algorithm was evaluated by comparing five reconstruction indicators: average running time (TIME), reconstruction error (ERROR), structural similarity (SSIM), peak signal-to-noise ratio (PSNR) and average number of iterations (Mean Iterations).

[0165] Table 1 is a comparison table of the four reconstruction indicators of reconstruction error (ERROR), average running time (TIME), structural similarity (SSIM) and peak signal-to-noise ratio (PSNR) for all images;

[0166] Table 1

[0167]

[0168] Table 2 is a comparison table of the mean iterations of all images;

[0169] Table 2

[0170]

[0171] Figure 4 This is a comparison chart of the reconstruction performance of 20 randomly selected sample images.

[0172] Table 1 summarizes the four reconstruction performance indicators of ERROR, TIME, PSNR, and SSIM on the three data sets. Each value in the table represents the average value of 50 images in each data set. It can be seen from Table 1 that the GAMP-FRSBL algorithm (the present invention) has better reconstruction performance. Its average error is the lowest on the three data sets, PSNR and SSIM are the highest, and reconstruction can be completed in a shorter time.

[0173] Table 2 lists the average number of iterations of several algorithms with better reconstruction performance on all images. It can be seen that compared with VGAMP-SBL and GGAMP-SBL, the GAMP-FRSBL algorithm can complete iterative convergence faster and significantly reduce the number of iterations and reconstruction time.

[0174] Figure 4 The reconstruction performance comparison chart of 20 randomly selected example images is listed. Figure 4 (a)-4(d) shows more intuitively that the GAMP-FRSBL algorithm has better reconstruction performance than other algorithms. Figure 4 (e) and 4(f), it can be seen that compared with the VGAMP-SBL and GGAMP-SBL algorithms, the GAMP-FRSBL algorithm converges faster, has a faster reconstruction time and better sparsity.

[0175] The results of the above embodiments confirm the rationality and effectiveness of the convex optimization method and damped GAMP algorithm used in steps S8-S10 of the present invention. The GAMP algorithm avoids matrix inversion, saves running time, and uses the damping coefficient to improve the convergence of the algorithm. The convex optimization method is used to iteratively update the hyperparameters, which improves the robustness of the model and the reconstruction accuracy. Its efficiency and performance are self-evident.

[0176] In summary, the present invention discloses a fast and robust sparse Bayesian image reconstruction method based on generalized approximate message passing. The present invention selects a public image data set, sets the pixel size of the image to 256×256 and processes it into a .mat format file to save the feature data of the image. Wavelet transform is added to each image for sparse preprocessing, and a multi-scale compressed sensing (CS) scheme is used for image reconstruction. The present invention first introduces the GAMP algorithm on the basis of SBL to replace the E step of the classic EM algorithm, avoids the matrix inversion problem, and greatly reduces the computational complexity of the model; and introduces a damping factor to enhance the convergence of the model. Finally, combined with the convex optimization strategy, the hyperparameters are iteratively updated by the block coordinate descent method (BCD) to improve the reconstruction accuracy of the model. The present invention improves the reconstruction efficiency and reconstruction accuracy by using the GAMP algorithm and the block coordinate descent method.

[0177] The present invention provides a fast and robust sparse Bayesian image reconstruction method based on generalized approximate message passing. There are many methods and ways to implement the technical solution. The above is only a preferred implementation method of the invention. All components not specified in this embodiment can be implemented using existing technologies.

[0178] The embodiments described above are merely descriptions of preferred implementations of the present invention and are not intended to limit the scope of the present invention. Without departing from the design spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by ordinary technicians in this field should all fall within the protection scope of the present invention.

Claims

1. A fast and robust sparse Bayesian image reconstruction method based on generalized approximate message passing, characterized by The steps include: S1, select image dataset; S2, set the pixel size of the image; S3, performing sparse preprocessing on the image; S4, establish a sparse Bayesian model; S5, calculating the likelihood function, the prior distribution of the sparse coefficient x, and the marginal likelihood function; S6, calculating the posterior density of the sparse coefficients by Bayes' theorem according to the likelihood function, prior distribution and marginal likelihood function calculated in step S5; S7, calculating the cost function; S8, iteratively update the hyperparameters using block coordinate descent method; S9, introduce the generalized approximate message passing algorithm to update the mean and variance; S10, introduce damping factor to enhance convergence; S11, judging whether the iteration termination condition is met, and if so, outputting the estimation of the sparse coefficient; If not satisfied, execute step S8; S12. Perform inverse wavelet transform on the output sparse coefficients to obtain a reconstructed image.

2. The fast robust sparse Bayesian image reconstruction method based on generalized approximate message passing as claimed in claim 1, characterized in that: In step S2, the pixel size of all images is set to 256×256 and processed into .mat format files to save the feature data of the images; and / or, in step S3, a "symmlet,8" wavelet transform with a coarsest scale of 3 and a finest scale of 7 is added to each image for sparse preprocessing.

3. The fast robust sparse Bayesian image reconstruction method based on generalized approximate message passing according to claim 1 or 2, characterized in that: In step S4, a sparse Bayesian model is established: y=Φx+e Where y∈R m is the observed data, x∈R n is the sparse coefficient, Φ∈R m×n is a known dictionary matrix, e∈R m is the noise vector; according to Bayes’ theorem, the estimate of x is obtained by maximizing the posterior distribution p(x|y; α, λ).

4. The fast robust sparse Bayesian image reconstruction method based on generalized approximate message passing as claimed in claim 3, characterized in that: In step S5, the likelihood function p(y|x; λ) is calculated according to the posterior distribution expression in step S4: Assume that the noise e has a mean of 0 and a variance of λ -1 Gaussian distribution of I, where I represents the identity matrix, i.e.: The hyperparameter λ follows the Gamma prior: Among them, shape parameter c>0, rate parameter d>0; is the Gamma function; The Gaussian likelihood function expression of the observed data y is obtained as: And / or, according to the posterior distribution expression in step S4, calculate the prior distribution p(x; α): in, is a diagonal matrix, α i is control x i The hyperparameter of sparsity, i=1,2,…,n; the hyperparameter α i Obey the Gamma prior distribution: Among them, shape parameter a>0, rate parameter b>0; is the Gamma function; And / or, according to the posterior distribution expression in step S4, calculate the marginal likelihood function p(y; α, λ): According to the properties of Gaussian function The marginal likelihood function p(y; α, λ) is approximately 0, and the covariance matrix is Gaussian distribution of; where the matrix F is a matrix, g is a vector, and C is a constant.

5. The fast robust sparse Bayesian image reconstruction method based on generalized approximate message passing as claimed in claim 4, characterized in that: In step S6, the posterior density p(x|y; α, λ) of x is calculated to obey the Gaussian distribution, and the posterior mean and the posterior covariance It is expressed as:

6. The fast robust sparse Bayesian image reconstruction method based on generalized approximate message passing as claimed in claim 5, characterized in that: step In S7, the cost function is calculated according to the posterior density obtained in step S8, as follows: λ and Λ α are all hyperparameters to be estimated, and α and λ are updated by type II maximum likelihood estimation by maximizing the joint probability density function p(y, α, λ): Maximizing the joint probability density function is equivalent to minimizing -2ln(α,λ); let σ 2 =λ -1 , γ i =α i -1 , Λ γ =diag(γ i )=Λ α ; Ignore the constant term and get the cost function L that combines x and hyperparameters x (x,γ,σ 2 )for:

7. The fast robust sparse Bayesian image reconstruction method based on generalized approximate message passing as claimed in claim 6, characterized in that: In step S8, the cost function L is minimized by iteration. x (x,γ,σ 2 ), and use block coordinate descent to adjust the hyperparameters γ and σ 2 To update:

8. The fast robust sparse Bayesian image reconstruction method based on generalized approximate message passing as claimed in claim 7, characterized in that: Step S9 In this paper, a generalized approximate message passing algorithm is introduced to update the mean and variance; For the prior p(x; α) of x and the likelihood function p(y|x; λ), the generalized approximate message passing algorithm uses Taylor expansion and quadratic approximation to approximate the maximum a posteriori estimate of x: Among them, q represents the system input vector and y represents the system output vector; therefore, the input function and output function are: Among them, the variable and Respectively Approximate estimates and noise-free models of The approximate values ​​of are and Let S = |Φ| 2 , S represents the square of the component amplitude of the dictionary matrix Φ; is the mean vector, Represents the covariance matrix The vector consisting of the diagonal elements of ; and represents the mean and variance of the intermediate variable s; therefore, the update of the input function and the output function is: The iterative update of the mean and variance is expressed as:

9. The fast robust sparse Bayesian image reconstruction method based on generalized approximate message passing as claimed in claim 8, characterized in that: In step S10, a damping factor is introduced to enhance the convergence of the model. The damping coefficient is calculated as follows: Among them, A0, B0 and C0 are constant terms; the damping factor θ s and θ x is the root of the above formula, θ s For control θ x For control 10. The fast robust sparse Bayesian image reconstruction method based on generalized approximate message passing as claimed in claim 9, characterized in that: The iteration termination condition of step S11 is: Determine whether the iterative convergence condition is met. If so, output the estimate of x: And / or, in step S12, the output sparse coefficient x is subjected to inverse wavelet transform to obtain a reconstructed image, and the calculation formula of the reconstruction error is: in, represents the reconstructed image, and P represents the real image.

Images