Disclosure of Invention
In view of the above, the present invention aims to provide a compressed sensing reconstruction algorithm (Hybrid Weighted Total Variation-non-local Rank, HWTV-NLR) based on a hybrid weighted total variation and a non-local low Rank, which mainly uses a local smooth sparse prior and a non-local self-similarity prior of an image, wherein the non-local similarity prior is used to preserve global information of the image, and the local smooth sparse prior suppresses noise while preserving detailed information of the image, so that a better effect can be achieved by effectively adjusting weights of the two items. In order to improve the quality of a reconstructed image, the invention improves the traditional total variation, namely, before reconstruction, an image high-frequency region (an edge, a noise and other regions) and a low-frequency region (a flat region) are detected through an edge detection operator, a second-order total variation is used for the low-frequency region of the image, and a first-order weighted total variation is used for the high-frequency region of the image.
In order to achieve the above purpose, the present invention provides the following technical solutions:
the image compressed sensing reconstruction method based on mixed weighted total variation and non-local low rank comprises the following specific processes:
s1: firstly, performing CS processing on an original image, and performing two-dimensional CS observation on the original image to obtain an observation value y; then decomposing the observed image by using an edge detector with an average curvature M structure to obtain a high-frequency part and a low-frequency part; the high frequency part mainly comprises edge details and noise of the image, and the low frequency part mainly comprises a flat area of the image; performing second-order total variation processing on the low-frequency part of the image, performing first-order weighted total variation processing on the high-frequency part of the image, and constructing a weight coefficient by a differential curvature C so as to obtain a mixed weighted total variation model;
s2: extracting sample blocks from the image, finding a certain number of sample similar blocks in a given area, and aggregating the image sample blocks and the sample similar blocks into a data matrix according to columns, wherein the data matrix has the property of low rank, so that a non-local low rank model is obtained;
s3: and combining the mixed weighted total variation model and the non-local low-rank model to obtain a final objective function, and carrying out optimization solution on the objective function by an alternate direction multiplier method.
Further, in the step S1, an edge detector is constructed by using an average curvature M based on a second order gradient, the image is adaptively decomposed into a high frequency part and a low frequency part, the low frequency part is processed by using a traditional second order total variation process, and the high frequency part is processed by using a first order weighted total variation process, wherein a weight coefficient is constructed by using a differential curvature edge detection operator based on the second order gradient, so that the edge detail of the image is protected, the robustness of the algorithm to noise is effectively improved, and finally, a new total variation model based on a total variation model, namely a mixed weighted total variation model is obtained, namely the expression of the mixed weighted total variation model is as follows:
wherein omegai As the weight coefficient, alphai For the edge detection operator, u represents the original reconstructed image, |Du|i Is first-order total variation, |D2 u|i Is a second-order total variation.
Further, in the step S1, the average curvature M is:
wherein,,ux and uy Respectively representing the first order of the current pixel point i of the original image along the horizontal direction and the vertical directionGradient uxx 、uxy 、uyx And uyy Representing the second order gradient of the current pixel point of the original image along different directions.
Further, in the step S1, the edge detector adaptively determines the high frequency and low frequency portions of the image by setting the threshold σ, when Mi At > sigma, alphai 1, wherein only a first-order total variation model acts on the high-frequency part; when Mi When less than or equal to sigma, alpha isi At 0, only the second-order total variation model acts on the smoothing portion.
Further, in the step S1, a weight coefficient ωi The method comprises the following steps:
wherein θ is a constant, Ci ＝||(uηη )i |-|(uξξ )i ||。
Further, in the step S3, the local smooth sparsity and the non-local self-similarity of the image are simultaneously depicted in the same frame, that is, the hybrid weighted total variation model is combined with the non-local low-rank model, and the combined model is used as a regular term constraint condition to construct an optimization model, that is, the objective function is expressed as:
the first term is a data fidelity term, the second term is a mixed weighted total variation regularization term, and the third term is a non-local low-rank regularization term; lambda (lambda)1 And lambda (lambda)2 As a regular term parameter, ωi R is the weight of each pixel pointi u represents a sample block ui And a matrix of similar image blocks, Li For a low rank matrix, Φ represents the observation matrix, y represents the observation value, η is the balance coefficient, and rank (·) represents the rank of the matrix.
The invention has the beneficial effects that: the invention considers the local smoothness and the non-local self-similarity of the image, improves the traditional mixed total variation model, constructs a weight coefficient by using a differential curvature by weighting the first-order total variation model, constructs a new edge detection operator by using an average curvature, and adaptively selects the total variation model by a threshold strategy. The mixed weighted total variation model combines the advantages of the first-order total variation model and the second-order total variation model, and effectively avoids the step effect while protecting the image edge. And redundant information is effectively reduced while the non-local low-rank model protects the global structural information of the image. In addition, the improved mixed weighted total variation model and the non-local low-rank model are taken as constraint construction prior regularization terms, and the soft threshold function and the smooth non-convex function are respectively adopted to solve the total variation and low-rank optimization problem, so that the structural characteristics of the image are fully utilized, and the adaptability and the noise immunity of the algorithm are effectively improved. Finally, the algorithm separates the target model into a plurality of sub-problems through an alternate direction multiplier method iteration strategy, and adopts the most efficient method to solve each sub-problem. Compared with the currently mainstream compressed sensing reconstruction algorithm, the algorithm has higher reconstruction quality, especially in the case of lower sampling rate. The method comprises the following steps:
1) The invention uses the mixed total variation to process the high frequency and low frequency of the image separately, processes the low frequency part of the image by the second-order total variation, suppresses the ladder effect, processes the high frequency part of the image by the first-order weighted total variation, suppresses the noise, and simultaneously keeps the edge details of the image, etc.
2) The invention introduces the mixed weighted total variation model into the non-local low-rank reconstruction algorithm, namely, the non-local self-similarity and the local smooth sparsity of the natural image are simultaneously and effectively depicted in a unified frame, and the adaptability and the reconstruction performance of the algorithm are enhanced.
Additional advantages, objects, and features of the invention will be set forth in part in the description which follows and in part will become apparent to those having ordinary skill in the art upon examination of the following or may be learned from practice of the invention. The objects and other advantages of the invention may be realized and obtained by means of the instrumentalities and combinations particularly pointed out in the specification.
Detailed Description
Other advantages and effects of the present invention will become apparent to those skilled in the art from the following disclosure, which describes the embodiments of the present invention with reference to specific examples. The invention may be practiced or carried out in other embodiments that depart from the specific details, and the details of the present description may be modified or varied from the spirit and scope of the present invention. It should be noted that the illustrations provided in the following embodiments merely illustrate the basic idea of the present invention by way of illustration, and the following embodiments and features in the embodiments may be combined with each other without conflict.
Referring to fig. 1 to 5, fig. 1 is a compressed sensing reconstruction algorithm based on hybrid weighted total variation and non-local low rank according to the present invention, and the specific processing flow is as follows:
(1) Firstly, an original image is subjected to two-dimensional CS observation to obtain an observation value y, wherein an observation matrix phi is a partial Fourier matrix, the observation value is subjected to inverse Fourier transform to obtain a pre-reconstructed image u, and each regularization parameter is initialized.
(2) For the non-local low-rank problem, firstly, finding out similar blocks of a sample image block according to a block matching method, wherein the used similar block matching method is based on Euclidean distance, then, aggregating the sample image block and similar blocks thereof into a matrix according to column arrangement, wherein the matrix has low-rank attribute, which is called a non-local low-rank matrix, and finally, converting the non-local low-rank problem into a minimum problem for solving rank, wherein the solving of the minimum problem for rank is approximated by using a non-convex proxy function log (°).
(3) For the total variation problem, a new edge detection operator is firstly constructed by utilizing an average curvature M based on a second derivative, an image is divided into a high-frequency part (such as an edge and noise) and a low-frequency part (a smooth area), and for the low-frequency area, second-order total variation processing is used; for the high frequency region, a first-order weighted total variation process is used, in which the first-order weighting coefficient is constructed by a differential curvature C, and the weight is smaller in the edge region and larger in the noise region. And finally, converting the total variation problem into a minimization optimization problem, and carrying out gradient solving of the high-frequency component and the low-frequency component by utilizing a soft threshold function.
(4) And combining the mixed weighted total variation and the non-local low-rank construction regularization term constraint to obtain a reconstructed objective function, decomposing the objective function into a plurality of sub-problems by using an Alternating Direction Multiplier Method (ADMM), and solving each sub-problem efficiently.
(5) And finally, carrying out experimental simulation by utilizing MATLAB, and carrying out comparative analysis by using the MATLAB and a mainstream image compressed sensing reconstruction algorithm to verify the effectiveness and feasibility of the algorithm.
1. Mixed weighted total variation model
In order to improve the adaptability and noise immunity of the traditional mixed total variation, the embodiment provides a new self-adaptive total variation strategy and a new first-order TV weighting model.
(1) Mixed total variation model
The traditional hybrid total variation model is generally expressed as:
wherein,,
wherein alpha isi For weighting coefficients, the edge detector selection is particularly important because the balanced hybrid total variation model is typically constructed using edge detection operators. The formula (2) and the formula (3) are a first-order total variation model and a second-order total variation model respectively, wherein (D)x u)i 、(Dy u)i The first-order difference operators of the current pixel point along the horizontal direction and the vertical direction are respectively adopted; (D)xx u)i 、(Dxy u)i 、(Dyx u)i 、(Dyy u)i And the second-order difference operators of the current pixel point along different directions are respectively adopted. Because of the complex image structure, the traditional edge detection operator based on the image gradient cannot effectively distinguish the edge, noise and smooth area of the image, so that the two regularization terms cannot be effectively balanced; and the first-order total variation model penalizes all gradients of the image uniformly, so that the edge of the image is protected and the denoising effect is poor in the reconstruction process. Aiming at the problems, a new edge detection operator is provided, a threshold value is set, a full-variation model is adaptively selected, and weights are set for the first-order full-variation model, so that the adaptability and the robustness of an algorithm are improved.
The present embodiment constructs a new edge detection operator using an average curvature M based on the second derivative, which can effectively distinguish between a high frequency part (such as an edge, noise) and a low frequency part (flat area) of an image, and the definition of the average curvature is as follows:
wherein,,ux and uy Represents the first-order gradient of the current pixel point i of the original image along the horizontal direction and the vertical direction, uxx 、uxy 、uyx And uyy Representing the second order gradient of the current pixel point of the original image along different directions. The first-order full-variation model can effectively protect the edge structure of the image, and the second-order full-variation model can effectively remove the step effect in the image smoothing area, so that the embodiment of the invention provides a self-adaptive full-variation strategy. Differentiating edges, noise mixed regions and smooth regions by setting a threshold sigma, when Mi At > sigma, alphai 1, only a first-order total variation model acts on an edge area and noise at the moment, and noise is removed while the edge of an image is protected; when Mi When less than or equal to sigma, alpha isi And 0, only the second-order total variation model acts on the smooth area, so that the step effect can be effectively avoided. The self-adaptive strategy fully utilizes the characteristics of the first-order total variation and the second-order total variation to achieve a better reconstruction effect.
(2) Hybrid weighted total variation
Aiming at the problem that all gradients of the image are punished to be consistent by the conventional unified full-variation model, the embodiment of the invention provides a novel mixed weighted full-variation model which can be expressed as:
wherein omegai As a weight coefficient, the traditional weighted total variation model only uses imagesThe weighting strategy provided by the invention utilizes a differential curvature C pixel edge detector based on a second-order gradient to construct a weighting coefficient, and can effectively distinguish an image edge region from a noise point, and the weighting coefficient at a pixel point i can be expressed as:
wherein,,
Ci ＝||(uηη )i |-|(uξξ )i || (7)
θ is a constant, at the edge, | (u)ηη )i Large, | (u)ξξ )i The I is very small, so C is very large, and omega is very small; at the noise point, | (u)ηη )i Sum (u)ξξ )i I are all large and almost equal, so C is small and ω is large. Therefore, the edge and the noise point of the image can be effectively distinguished according to the differential curvature C, and the first-order total variation model has smaller weight in the edge area, so that the edge of the image can be effectively protected; the first-order total variation model has larger weight at noise points, so that noise can be removed better, and the robustness of the model is improved effectively by the novel weighting strategy.
2. Non-local low rank model
Many repeated similar structures exist in a natural image, and for any sample image block, similar blocks similar to the image can be found in the image, and the property is called non-local similarity of the image. The construction of the non-local low rank matrix is as follows:
assume that the sample block size at pixel i position isDenoted as ui ∈Rn ，ui Looking up m image blocks similar to the local window (20×20):
wherein G isi Represents the positions in the image of all similar blocks of sample block i, d () represents the euclidean distance, T represents the similarity threshold,representing the j-th similar block of sample block i.
Finding ui After each similar block of (2), each similar block is sequentially expanded in sequence from high to low in similarity, so that a data matrix X is obtainedi ＝[ui,0 ui,1 …ui,m-1 ]，Xi ∈Rn×m . Since these image blocks have similar structural information, a data matrix X is formedi With low rank properties. In practice, Xi May be affected by noise, one solution is to matrix the data Xi Modeled as Xi ＝Li +Si Wherein L isi Is a low rank matrix, Si Is a gaussian noise matrix, and therefore, the NLR-CS algorithm model can be expressed as:
wherein,,representing the variance of additive gaussian noise. Converting it into unconstrained form is:
wherein R isi u＝[Ri1 u,Ri2 u,…Rih u]Representing a sample block ui And a matrix of similar image blocks.
3. Combined model
The objective function of the above-described mixed weighted total variation and non-local low rank CS image reconstruction algorithm is as follows:
wherein lambda is1 And lambda (lambda)2 As a regular term parameter, ωi The weight of each pixel point. The first term is a data fidelity term, the second term is a mixed weighted total variation regularization term, and the third term is a non-local low-rank regularization term.
Because the formula (11) contains low-rank matrix minimization and different gradient optimization problems, direct solution is very complex, in order to simplify the solution of an algorithm model, the embodiment of the invention adopts an ADMM method framework to carry out optimization solution on the proposed objective function, and the ADMM method has the advantages that when the model meets certain convexity conditions, the model can be converged to the optimal solution of the original model without meeting the condition that parameters tend to infinity, and the original problem is solved separately and independently, so that the solution of the model becomes easier. The method comprises the following specific steps:
by introducing auxiliary variables x, z1 And z2 The unconstrained minimization problem of equation (11) is rewritten as a constraint problem as follows:
s.t.u＝x,Du＝z1 ,D2 u＝z2
the corresponding augmented lagrangian function is:
wherein mu1 、μ2 Sum mu3 Penalty parameters for controlling linear constraints; a. b and c are lagrange multipliers; the ADMM method can be regarded as an original dual method of solving the saddle point problem (augmented lagrangian), soThe saddle point of equation (13) is the optimal solution of the original problem (11). Next, we decompose equation (13) into a number of simpler minimisation sub-problems a-E, resulting in a corresponding solution.
A、z1 -sub-problem
Fix u, Li X, and z2 Then z1 The optimization problem of (2) can be expressed as:
the sub-problem can be solved by a soft threshold operator, and can be obtained:
wherein sign (·) is a sign function, and when·is equal to or greater than 0, sign (·) =1; when · is less than 0, sign (·) = -1.
B、z2 -sub-problem
The same as the A sub-problem, fix u, Li X, and z1 Then z2 The optimization problem of (2) can be expressed as:
the sub-problem can also be solved by a soft threshold operator, and can be obtained:
wherein I is an identity matrix.
C、Li -sub-problem
Fix u, x, z1 And z2 L is theni The optimization problem of (2) can be expressed as:
since the low rank solution of the matrix is an NP-hard problem, it cannot be directly solved, based on which Dong et al uses a smooth and non-convex function L (Li Epsilon) to approximate rank (Li ) And effectively proves that the low-rank approximation effect of logdet () is better than the core norm I I.I.I.)* (sum of singular values). L (L)i Epsilon) can be defined as:
L(Li ,ε)＝logdet((Li LiΤ )1/2 +γI) (19)
wherein, gamma is a positive constant, I is an identity matrix, and the iterative process is as follows:
wherein,,is Xi Is a singular value decomposition of τ=η/2λ2 ，(x)+ ＝max(x,0)，ξj(k) ＝1/(σj(k) +ε)，σj Is Li The j-th singular value.
D. x-child problem
Fix u, Li 、z1 And z2 The optimization problem for x can be expressed as:
equation (21) is essentially a strict convex quadratic function minimization problem, so let xk+1 The derivative of (2) is zero to obtain its closed solution as shown below
Wherein I is an identity matrix.For a diagonal matrix, each element in the diagonal matrix corresponds to an image pixel location, and the value of the element is the number of overlapping patches covering the pixel location. />Representing the average result of the image blocks, all collected similar blocks are averaged for each sample block.
E. u-child problem
Fix Li 、z1 、z2 And x, then the optimization problem for u can be expressed as:
equation (23) is also a strict convex quadratic optimization problem containing a closed-form solution, and the derivative is zero, so that
To solve equation (24) more efficiently, we solve using the Preconditioned Conjugate Gradient (PCG) method.
After solving the above five sub-problems, the lagrangian multiplier is updated according to equation (25), as follows:
then, the next iteration is performed until the iteration termination condition is satisfied. Output u of last iterationk+1 Reconstructed image that is original image
In order to verify the effectiveness and feasibility of the invention, we performed experimental simulations on the inventive examples. Six images were selected for simulation, barbara, boats, cameraman, lena, monarch, parrots each with a resolution of 256×256. In order to more objectively verify the invention, the most mainstream CS image reconstruction algorithm at present is selected for experimental comparison, and the experimental comparison is respectively as follows: NLR-CS algorithm, TVNLR algorithm. The performance evaluation index is PSNR (peak signal to noise ratio, dB). All simulation experiments were performed on Matlab2016 a.
In the comparative experiment, the settings of the relevant parameters were as follows: m is initialized to be an identity matrix; k, T, K0 18,15 and 5, respectively; the lagrangian multipliers a, b, c are set according to sampling rates, different sampling rates possibly corresponding to different values; selecting the most proper iteration times to be 270 times according to experimental simulation; the number of sample-like blocks is 45; the size of the similar block matching window is set to 20×20, and the window is centered on the current pixel point; the sample block has a size of 6×6.
Fig. 2 shows analysis of stability of the present invention at sampling rates of 0.05 (fig. 2 (a)) and 0.20 (fig. 2 (b)), and it can be seen from the graph that, under different sampling rates, the PSNR curve monotonically increases with increasing iteration number, and when the iteration number reaches a certain number, the PSNR curve finally becomes stable, and the experiment fully proves the stability of the algorithm of the present invention.
Fig. 3 and fig. 4 are subjective visual contrast diagrams of simulation results of Barbara, camaraman images at a sampling rate of 0.05, respectively. The image reconstruction method has the advantages that the algorithm provided by the invention obtains the maximum PSNR value under the same condition, the image reconstruction effect is the best, and the effectiveness of the algorithm provided by the invention is verified; in order to verify the robustness of the proposed algorithm to noise, in the simulation process, different levels of noise are added to the observed value to carry out experimental simulation, and the reconstruction effects of the proposed algorithm and the comparison algorithm under the noise condition are compared. Taking the sample rate of 0.05 as an example, the noise is gaussian white noise, and fig. 5 (a) and fig. 5 (b) are reconstruction performance curves of the images barbera and Cameraman at different noise levels, respectively. As can be seen from fig. 5, the PSNR values of the proposed algorithm are higher than those of the comparative algorithm in the case of being interfered by noise of different degrees, which verifies that the algorithm has better robustness to noise.
Finally, it is noted that the above embodiments are only for illustrating the technical solution of the present invention and not for limiting the same, and although the present invention has been described in detail with reference to the preferred embodiments, it should be understood by those skilled in the art that modifications and equivalents may be made thereto without departing from the spirit and scope of the present invention, which is intended to be covered by the claims of the present invention.