CN107784278B - Sparse image reconstruction accuracy reduction complexity method is improved with structuring is prior-constrained - Google Patents
Sparse image reconstruction accuracy reduction complexity method is improved with structuring is prior-constrained Download PDFInfo
- Publication number
- CN107784278B CN107784278B CN201710969340.9A CN201710969340A CN107784278B CN 107784278 B CN107784278 B CN 107784278B CN 201710969340 A CN201710969340 A CN 201710969340A CN 107784278 B CN107784278 B CN 107784278B
- Authority
- CN
- China
- Prior art keywords
- matrix
- value
- wavelet
- parameter
- coefficient
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/22—Source localisation; Inverse modelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
- G06F18/2134—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on separation criteria, e.g. independent component analysis
- G06F18/21345—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on separation criteria, e.g. independent component analysis enforcing sparsity or involving a domain transformation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/02—Preprocessing
- G06F2218/04—Denoising
- G06F2218/06—Denoising by applying a scale-space analysis, e.g. using wavelet analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/08—Feature extraction
- G06F2218/10—Feature extraction by analysing the shape of a waveform, e.g. extracting parameters relating to peaks
Abstract
Sparse image reconstruction accuracy reduction complexity method is improved with structuring is prior-constrained, belong to Remote Sensing Image Compression technical field, it is characterized in that, it is a kind of method for introducing image information on the basis of more observation vector models to improve precision reduction complexity, successively the following steps are included: input observing matrix, the matrix of more observation vector compositions, the wavelet transformation number of plies, every layer of wavelet coefficient number；Initialize system model parameter, image wavelet coefficient, maximum number of iterations, precision property index；Image wavelet coefficient is reconstructed with the gibbs sampler based on Markov chain monte carlo method, then wavelet inverse transformation is carried out to matrix of wavelet coefficients after reconstruct, to obtain reconstructed image.The present invention can apply required precision property index to practical reconstructed image, complexity be reduced using more observation vector models, while adjusting reconstruction accuracy further through the method for the adjustment wavelet transformation number of plies, so that the quality of reconstructed image meets performance indicator requirement.
Description
Technical field
The invention belongs to field of signal processing, can be applied to the sparse reconstruct of image information, in particular to are seen based on more
Survey the low complex degree structuring reconstructing method of vector model.
Background technique
According to compressed sensing (Compressive Sensing, CS) theory, if signal meets centainly in a certain transform domain
Sparsity can then sample signal using the rate far below nyquist sampling rate, and realize to the accurate of signal
Reconstruct.Currently, compressed sensing technology has application, such as Sub-nyquist sampling system, compression imaging in many engineering fields
System, compression sensing network etc..Particularly, compressed sensing technology can also be applied on sparse image reconstruction, and
Through being quite widely applied in key areas such as medical imagings.
However, existing compressed sensing algorithm is mostly based on single observation vector (SMV) when solving the problems, such as image reconstruction
Entire sparse image is considered as one-dimensional vector and handled by model；The problems of the processing mode is: when our needs
When high-definition picture is reconstructed, the one-dimensional vector dimension being drawn by large-scale image is quite high, therefore relates in algorithm
And the observing matrix and covariance matrix arrived is in large scale, so that sizable storage pressure can be brought to existing computing platform
With calculating pressure.For example, if opening the image of 1024 × 1024 dimensions using the compressed sensing algorithm process one based on SMV,
Dimension will be more than 1.04 × 10 after the image array being related in algorithm transforms into a column vector6Dimension, stores such scale
Covariance matrix need be more than 64GB computing platform memory；Even if computing platform memory is sufficiently large, matrix storage can satisfy
Demand simultaneously makes these methods can be with successful operation, but the computing resource for the calculating time and occupancy thus expended is also to be difficult to hold
It receives.Therefore, in high-definition picture reconstruction, do not have just with the compressed sensing algorithm based on SMV model existing
Sincere justice.
For this purpose, the present invention fully considers the characteristics of image sparse reconstruct, problem is characterized again using MMV model.
MMV model handles each column vector of original signal as independent subsignal, and in array signal processing, cognition
The fields such as radio have been achieved for corresponding research achievement, but are not applied to the sparse reconstruction of image.Therefore, such as
More observation vector models are introduced large-scale image reconstruction to fruit by us, and each column vector of image array are regarded as
It is an observation vector in MMV model, then such model characteristic manner does not further relate to convert image to high-dimensional single
The operation of vector, corresponding problem scale will be reduced significantly, so that being solved using the principle of compressed sensing big
The reconstruct of scale image becomes the engineering problem that can be realized in most of computing platforms.But if by existing solution
The MMV restructing algorithm of other field problem directly applies to when solving the problems, such as image reconstruction, then having ignored image information is had
Architectural characteristic, to cause the serious loss and the deterioration of reconstruction accuracy of information.
It therefore, is the reconstruction accuracy deterioration problem being likely to occur when further overcoming directly using MMV model, institute of the present invention
Reconstructing method is stated on the basis redescribed using MMV model to problem, introduces in image information and potentially ties
Structureization is prior-constrained, be embodied as of both structuring prior information: on the one hand, image it is dilute on wavelet transformed domain
Sparse coefficient has the probability dependency of quaternary tree shape structure, on the other hand, between the column vector of image array and dependent, and
It is that there is certain spatial coherence.
In conclusion the present invention is by introducing high-definition picture reconstruction for MMV model, can be realized makes its storage
Complexity and computation complexity decline to a great extent；Further, prior-constrained by introducing the structuring of image information, it ensure that
The precision of image reconstruction under MMV model.
Summary of the invention
(1) technical problems to be solved
The storage complexity faced present invention aim to address the sparse reconstruct of high-definition picture and time are complicated
Spend the limited problem of higher and reconstruction accuracy.
(2) technical solution
In order to solve the above-mentioned technical problem, the present invention proposes a kind of structure based on more observation vectors characterization of low complex degree
Change reconstructing method.The purpose of the present invention is what is be achieved through the following technical solutions: input observing matrix, more observation vector compositions
Matrix, the wavelet transformation number of plies, every layer of wavelet coefficient number；Initialize system model parameter, image wavelet coefficient, greatest iteration time
Number, precision property index；Weight is carried out to image wavelet coefficient with the gibbs sampler based on Markov chain monte carlo method
Structure, then wavelet inverse transformation is carried out to matrix of wavelet coefficients after reconstruct, so that reconstructed image is obtained,
1, sparse image reconstruction accuracy reduction complexity method is improved with structuring is prior-constrained, which is characterized in that be one
Kind is the remote sensing images X of N × N=1024 × 1024 to resolution ratio0After being redescribed using more observation vector models, then benefit
Use X0Sparse coefficient on wavelet transformed domain has the probability dependency and more observation vector matrixes of quaternary tree shape structure
The structuring of spatial coherence these two aspects between column vector it is prior-constrained come reconstruct resolution ratio be N × N=1024 ×
1024 remote sensing images X*Method, successively realize according to the following steps in a computer:
Step (1) inputs following parameter to computer:
The remote sensing images matrix X of N × N=1024 × 10240, abbreviation matrix X0,
It is generated with the tool box MATLAB and is based on the matrix X0Quaternary tree shape structure 2-d discrete wavelet transformation of coefficient
Matrix Θ, resolution ratio are N × N=1024 × 1024, and the number of plies S of two-dimensional discrete wavelet conversion is based on the matrix X0Point
Resolution is N × N=1024 × 1024, the maximum value S of the number of pliesmax=8, level number s=1,2 ..., 8, quaternary tree shape structure
The number of each layer of wavelet coefficient is (NC) in wavelet coefficients, root node is in first layer S1, appoint and take S=4,
With the tool box MATLAB generate M × N random observation matrix, M≤N, observing matrix abbreviation matrix Φ, similarly hereinafter,
The matrix Y of more observation vectors composition based on the matrix Θ M × N generated, this is to be with the matrix Φ
For the matrix X under conditions of observing matrix0Observation, i.e. Y=Φ Θ,
The remote sensing images X obtained after input reconstruct*Two image property indexs: the setting value of Y-PSNR PSNR,
The setting value of mean square error MSE,
Step (2) initializes Modulus Model parameter, comprising:
1. the regularized covariance Matrix C of the wavelet coefficient of different level numbers1, C2..., C4, when initialization is all unit square
Battle array, abbreviation Cs={ C1, C2..., C4,
2. controlling each CsThe parameter lambda of value, enables λ be initialized as 2,
3. the precision parameter for the Gaussian Profile that observation noise is obeyed is the α reciprocal of the variance of observation noisen, at the beginning of enabling it
Beginning turns to αn=10,
4. the precision parameter α for the Gaussian Profile that the wavelet coefficient of different level numbers is obeyed1, α2..., α4, it is each layer small echo
The inverse of variance between coefficient column vector, when initialization, are equal to 1,
5. setting: control the probability that each layer wavelet coefficient takes zero:
For first layer S1: control matrix of wavelet coefficients Θ(1)In elementValue take zero probability to be initialized asAnd enable W1=0.3, correspondingly,The probability of negated zero is 1-W1, wherein symbol Θ(1)、In on be designated as the level number of each layer matrix of wavelet coefficients Θ, subscript (ij) is sequence of each element W in row, column
Number, similarly hereinafter,
For the second layer to the 4th layer of S2, S3, S4:
If matrix of wavelet coefficients Θ(s), the element of 2≤s≤4The coefficient value of its father's node under quad-tree structureThen controlValue takes zero probability2 s≤4 <, and initialize W0s=0.3, symbol W0sIn
Symbol " 0 " indicate father's nodeCoefficient value be zero, symbol " π " indicate father's node, similarly hereinafter,
If matrix of wavelet coefficients Θ(s), the element of 2≤s≤4The coefficient value of its father's node under quad-tree structureThen controlValue takes zero probability2 s≤4 <, and initialize W1s=0.7, symbol W1sIn
Symbol " 1 " indicate father's nodeCoefficient value be 1,
6. system hyper parameter a0,b0,c0,d0,e0,f0, the hyper parameter is that the parameter in a kind of pair of step (2) is obeyed
The parameter obtained after prior probability distribution parametrization, in which:
a0,b0,c0,d0In initialization all equal to 10-6, it is constant for controlling the parameter of probability-distribution function,
Symbol e0It include: (ec)1, and (ecz)s, 2≤s≤4:
(ec)1For controlling W1The probability density function obeyed, value are as follows:
(ec)1=0.9 × (NC)1,
(ecz)2, (ecz)3, (ecz)4, it is respectively used to W0 when 2≤s of control≤42, W03, W04The probability density obeyed point
Cloth function, value are as follows:
(ecz)2=0.1 × (NC)2, (ecz)3=0.1 × (NC)3, (ecz)4=0.1 × (NC)4, (eco)2, (eco)3,
(eco)4It is respectively used to W1 when 2≤s of control≤42, W13, W14The probability density function obeyed, value are as follows:
(eco)2=0.5 × (NC)2, (eco)3=0.5 × (NC)3, (eco)4=0.5 × (NC)4,
Symbol f0It include: (fc)1, and (fcz)s, 2≤s≤4:
(fc)1For controlling W1The probability density function obeyed, value are as follows:
(fc)1=0.1 × (NC)1,
(fcz)2, (fcz)3, (fcz)4It is respectively used to W0 when 2≤s of control≤42, W03, W04The probability density obeyed point
Cloth function, value are as follows:
(fcz)2=0.9 × (NC)2, (fcz)3=0.9 × (NC)3, (fcz)4=0.9 × (NC)4, (fco)2, (fco)3,
(fco)4It is respectively used to W1 when 2≤s of control≤42, W13, W14The probability density function obeyed, value are as follows:
(fco)2=0.5 × (NC)2, (fco)3=0.5 × (NC)3, (fco)4=0.5 × (NC)4,
7. the parameter t, T of initialization control iterative processmax, in which:
T is the serial number of the number of iterations, TmaxFor maximum number of iterations, thus the setting value of the serial number of the number of iterations be t=1,
2,...,Tmax}；
Step (3) successively calculate by following procedure the iterative process of matrix of wavelet coefficients Θ:
Step (3.1) successively carries out first time iteration according to the following steps, enables t=1,
Step (3.1.1) sets various system model parameters in constraint condition and step (2):
Constraint condition are as follows: the remote sensing images matrix X of N × N=1024 × 10240, the wavelet systems of N × N=1024 × 1024
Matrix number Θ is quad-tree structure, Smax=8, appoint and takes S=4,
When initial, the matrix of wavelet coefficients Θ of each layer is full null matrix,
The initial value according to system model parameter and wavelet coefficient is successively calculated as follows in step (3.1.2), calculates intermediate
Variable
Wherein:
||·||2Indicate the 2- norm of amount of orientation or matrix,
The element Θ arranged for the i-th row of Θ matrix, jthijThe precision parameter for the probability distribution obeyed,
The element Θ arranged for the i-th row of Θ matrix, jthijThe mean of a probability distribution parameter obeyed,
The element Θ arranged for the i-th row of Θ matrix, jthijThe adjustment parameter of zero probability is taken,
αsAnd αnSystem model parameter value when expression is initial,System model when expression is initial
The intermediate variable that parameter is calculated,
Φ-iIndicate to remove the matrix of other all column compositions except its i-th column in Φ,
Φ·iIndicate the vector of all elements composition of the i-th column of Φ matrix,Indicate vector Φ·iTransposition,
Θ(-i)jIt indicates in vector theta·jThe vector of the middle other all elements composition removed except its i-th of element,
Step (3.1.3), according to intermediate variable obtained in step (3.1.2)Calculate wavelet coefficient Θ institute
Posterior probability distribution p (the Θ of obedienceij):
Wherein: δ0The uni-impulse function at origin is represented,
Posterior probability distribution p (the Θ obeyed according to wavelet coefficient Θij), sampling obtains the i-th row of matrix Θ, jth column
Element ΘijValue to get arrive matrix Θ sampled value,
Step (3.1.4), according to the sampling of matrix Θ obtained in the initial value of system model parameter and step (3.1.3)
Value calculates system model parameter is obeyed in first time iteration posterior probability Density Distribution and Matrix CsFirst time change
For result:
Wherein:
The value set of footmark s is { 1 ..., 4 }, it is noted that only in the S of input >=2, and the number of plies s currently considered
When greater than 1, just there is W0sAnd W1s,
The Frobenius norm of representing matrix,
||·||0The 0- norm of representing matrix (or vector),
Θ(s)It is that s layers of wavelet conversion coefficient combines the matrix to be formed according to natural order,
RsIt is ΘsLine number,Indicate ΘsR row,
Θs,kIndicate s layers of k-th of wavelet coefficient, Θπ(s,k)It is Θs,kFather's node coefficient value,
P (x)=Gamma (a, b) indicates that stochastic variable x obeys the gamma that parameter is (a, b) and is distributed, and expanded form isWherein Γ () indicates gamma function, by taking a is independent variable as an example, the expanded form of gamma function
For
P (x)=Beta (a, b) indicates that stochastic variable x obeys the beta that parameter is (a, b) and is distributed, and expanded form is
Indicating indicative function, i.e., functional value takes 1 when the expression formula in bracket is true, otherwise functional value takes 0,
According to above system model parameter αn、αs、W1、W0sAnd W1sThe posterior probability Density Distribution obeyed, to system mould
Shape parameter is sampled, and updates corresponding system model parameter value, for iterative process next time,
Step (3.1.5), the sampled value of matrix Θ obtained in storing step (3.1.3), first time iterative process are formal
Terminate,
Step (3.2) carries out second to T according to the description of step (3.1.1)~step (3.1.5)maxSecondary iteration
Calculating,
Step (3.3), to the serial number of all the number of iterations, t={ 1 ..., Tmax, it calculates storage in each iteration and obtains
All matrix Θ arithmetic mean of instantaneous value, obtained average value is denoted as Θ*, to average value Θ*Carry out S layers of 2-d wavelet inversion
It changes, obtains X*,
Step (3.4), counts reconstructed image X*Relative to original image X0Error extension, PSNR and MSE, if this
Two indexes meet the performance indicator requirement of input, then algorithm stops, the wavelet coefficient number of plies otherwise inputted in set-up procedure (1)
S, and S:=S+1 is enabled, go to step (2) again, again operation initialization and iterative process, until reconstructed image is opposite
Meet performance indicator requirement in the error extension of original image.
Beneficial effect
The present invention gives a kind of sparse reconstructing methods of the high-definition picture of low complex degree.The reconstructing method is based on
MMV model is designed；To further increase reconstruction accuracy, the structuring that this method introduces image is prior-constrained.This method
Have the advantage that (1) makes the storage complexity of problem by O (N by introducing MMV model4) it is reduced to O (N2), it calculates complicated
Degree is by O (M2N2) it is reduced to O (MN), to provide feasible scheme for high-definition picture reconstruction；(2) pass through introducing
The constraint of structuring prior probability distribution, improves the accuracy of image reconstruction under MMV model；(3) it utilizes and is derived in this method
The Gibbs sampling method arrived, can conclude that obtain image coefficient and model parameter under the structuring MMV model it is theoretical most
Excellent solution.
Detailed description of the invention
Fig. 1 is the flow chart that the present invention realizes sparse image reconstruction.
Fig. 2 features the quaternary tree shape probability dependency structure of image multilayer wavelet transform, and wherein s is that multi-level Wavelet Transform becomes
The number of plies changed.
Fig. 3 illustrates the stratification Bayesian model based on structuring prior information that the present invention establishes.Wherein a0,b0,
c0,d0,e0,f0The hyper parameter of picture engraving coefficient and model parameter prior distribution, Θ, Y be respectively image wavelet coefficient and
Based on the sparse observing matrix that more observation vector models obtain, αs,Cs,ωsIt is the parameter for portraying the probability density characteristics of Θ, αn
It is the parameter for portraying the probability density characteristics of Y.
Fig. 4 is that image reconstructing method proposed by the present invention (cM-SSBL) and the existing sparse reconstruct based on MMV model are calculated
Contrast on effect of the method (OMP-MMV, M-FOCUSS and MSBL) in the sparse reconstruction of processing large-scale image.Horizontal axis is normalizing
Change sample rate M/N, the longitudinal axis is Relative reconstruction error；It indicates to use OMP-MMV algorithm,Using M-FOCUSS
Algorithm,It indicates to use MSBL algorithm,It indicates to use cM-SSBL algorithm.
Fig. 5 is that image reconstructing method proposed by the present invention (cM-SSBL) and the existing sparse reconstruct based on MMV model are calculated
Contrast on effect of the method (OMP-MMV, M-FOCUSS and MSBL) in the sparse reconstruction of processing large-scale image.Horizontal axis is normalizing
Change sample rate M/N, the longitudinal axis is Y-PSNR (PSNR)；It indicates to use OMP-MMV algorithm,Using M-
FOCUSS algorithm,It indicates to use MSBL algorithm,It indicates to use cM-SSBL algorithm.
Fig. 6 is to use resolution ratio for 1024 × 1024 standard testing image " cameraman ", is in normalization sample rate
In the case where 0.4, image is carried out to OMP-MMV, M-FOCUSS, MSBL and tetra- kinds of methods of cM-SSBL proposed by the invention
The visual effect display diagram of reconstitution experiments.Fig. 6 (a) is the reconstruction result of OMP-MMV algorithm and corresponding PSNR, Fig. 6 (b) are M-
The reconstruction result of FOCUSS algorithm and corresponding PSNR, Fig. 6 (c) are the reconstruction result and corresponding PSNR, Fig. 6 of MSBL algorithm
(d) be cM-SSBL algorithm proposed by the present invention reconstruction result and corresponding PSNR.
Fig. 7 is the reference curve during parameter preferably relevant to iteration control, and horizontal axis is Y-PSNR, the longitudinal axis
The number of iterations, in order to show the influence of the variation of the number of iterations to reconstruct error performance, Fig. 7 intercepted preferred process carry out to
Finally, Y-PSNR situation of change when the number of iterations is relatively fewer.
Specific embodiment
With reference to the accompanying drawings and examples, specific embodiments of the present invention will be described in further detail.Implement below
Example is not limited the scope of the invention for illustrating the present invention.
As shown in Figure 1, the purpose of the present invention is be achieved through the following technical solutions: input observing matrix, observe more to
Measure the matrix of composition, the wavelet transformation number of plies, every layer of wavelet coefficient number；Initialization system model parameter, image wavelet coefficient, most
Big the number of iterations, precision property index；With the gibbs sampler based on Markov chain monte carlo method to image wavelet system
Number is reconstructed, then carries out wavelet inverse transformation to matrix of wavelet coefficients after reconstruct, to obtain reconstructed image.The present invention can be to reality
Border reconstructed image applies required precision property index, reduces complexity using more observation vector models, while further through adaptive
The method adjustment the number of iterations answered, so that the precision of reconstructed image meets performance indicator requirement.
Now for carrying out sparse reconstruct to standard testing image X of the width dimension for N × N (N=1024),
Sparse observing matrix is determined first, and the observation of M × N (M=512) is generated using MATLAB generating random number tool box
Matrix Φ,
According to the basic principle that image sparse characterizes, for original remote sensing images under conditions of matrix Φ is characterization substrate
X0It is observed, obtains the matrix Y of M × N of observation vector composition,
The number of plies S=4 of two-dimensional discrete wavelet conversion is set, then the number of every layer of wavelet coefficient also determines therewith: in N=
In the case where 1024, S=4, the 0th layer of wavelet coefficient number is (NC)0=3 × 642=12288, the 1st layer of wavelet coefficient number be
(NC)1=3 × 1282=49152, the 2nd layer of wavelet coefficient number is (NC)2=3 × 2562=196608, the 3rd layer of wavelet coefficient
Number is (NC)3=3 × 5122=786432；
(this sentences MSE=16.3 setting performance indicator to setting reconstructed image performance indicator PSNR=36, and effect is of equal value
)；
Next, the sparse image reconstruction algorithm based on structuring prior information proposed according to the present invention, successively carries out
Following solution procedure:
The wavelet coefficient Θ for initializing N × N enables Θ be equal to full null matrix；
Initialize system model parameter, comprising:
The regularized covariance Matrix C of the wavelet coefficient of different level numbers0,C1,...,CS-1, unit matrix is all when initial,
Abbreviation CS,CS={ C0,C1,...,CS-1}；
Control CSThe parameter lambda of size, enables λ=2；
The precision parameter for the Gaussian Profile that observation noise is obeyed is the inverse of the variance of observation noise, is denoted as αn, αn=
10；
The precision parameter α for the Gaussian Profile that the wavelet coefficient of different level numbers is obeyed0,α1,...,αS-1, it is s layers respectively
The inverse of covariance between wavelet coefficient column vector, when initialization, are equal to 1；
W1Control first layer wavelet coefficient takes zero probability, and sequence s, s=2 ..., S, W0 are designated as under2...,
W0SAnd W12..., W1SThe probability that wavelet coefficient takes zero is controlled when being all s > 1,
For scale be N × N wavelet coefficient Θ the i-th row corresponding at different sequence s, jth column member
The setting value of element, i=1,2, i .N, j=1,2, j., N:
W1Θ when for first layer s=1ij|S=1WhenValue,
W0sFor s > 1 and Θπ(ij)When=0Value,Wherein,
Θπ(ij)It is the Θ in the case where Θ is quad-tree structureijFather's node coefficient value, similarly hereinafter, s value range be { 2 ..., S }, i=
1,2, i .N, j=1,2, j., N,
W1sFor s > 1 and Θ π(ij)When=1Value,S value range
For { 2 ..., S }, i=1,2, i .N, j=1,2, j., N,
Initialization system hyper parameter a0,b0,c0,d0,e0,f0, in which:
a0,b0,c0,d0In initialization all equal to 10-6,
Symbol e0It include: (ec)1, and (ecz)2..., (ecz)s(eco)2..., (eco)s:
(ec)1For controlling W1The probability density function obeyed, value are as follows:
(ec)1=0.9 × (NC)1,
(ecz)2..., (ecz)SFor controlling W02..., W0SThe probability density function obeyed, value are as follows:
(ecz)2=0.1 × (NC)2..., (ecz)S=0.1 × (NC)S, (eco)2..., (eco)SFor controlling
W12..., W1SThe probability density function obeyed, value are as follows:
(eco)2=0.5 × (NC)2..., (eco)S=0.5 × (NC)S,
Symbol f0It include: (fc)1, and (fcz)2..., (fcz)S(fco)2..., (fco)S:
(fc)1For controlling W1The probability density function obeyed, value are as follows:
(fc)1=0.1 × (NC)1,
(fcz)2..., (fcz)SFor controlling W02,...,W0S-1The probability density function obeyed, value is such as
Under:
(fcz)2=0.9 × (NC)2..., (fcz)S=0.9 × (NC)S,
(fco)2..., (fco)SFor controlling W12,...,W1S-1The probability density function obeyed, value is such as
Under:
(fco)2=0.5 × (NC)2..., (fco)s=0.5 × (NC)s；
The parameter t, T of initialization control iterative processmax, in which:
T is the serial number of the number of iterations,
Set maximum number of iterations Tmax=55.
Successively calculate by following procedure the iterative process of wavelet coefficient Θ,
Step (1) determines the initial value of system model parameter according to current wavelet coefficient number of plies S；
Step (2), TmaxFor the good constant of current setting, the periodicity at most iterated to calculate is represented,
On be designated as the number of iterations,
Step (3) successively realizes iterative calculation for the first time, serial number t=t as follows(1)=1,
Initial value of the step (3.1) using the initial value of each system model parameter and wavelet coefficient as first time iteration,
The initial value according to system model parameter and wavelet coefficient is successively calculated as follows in step (3.2), calculates intermediate become
Amount
Step (3.3) calculates the Posterior probability distribution p (Θ that wavelet coefficient Θ is obeyedij):
Step (3.4), the Posterior probability distribution p (Θ obeyed according to wavelet coefficient Θij), sampling obtains the of matrix Θ
I row, jth column element ΘijValue to get the sampled value for arriving matrix Θ, according to the initial value of system model parameter and matrix Θ
Sampled value calculates system model parameter is obeyed in first time iteration posterior probability Density Distribution and Matrix CsFirst
Secondary iteration result,
Step (4) carries out second to T according to the description of step (3.1)~step (3.4)maxThe calculating of secondary iteration,
Step (5), to all Θ(t),Calculate all Θ(t)Arithmetic mean of instantaneous value, what is obtained is flat
Mean value is denoted as Θ*, to average value Θ*The wavelet inverse transformation for carrying out S layers, obtains X*,
Step (6), counts reconstructed image X*PSNR, do not meet the performance indicator of PSNR=36 yet at this time, then by Tmax
Increase 10, at this time Tmax=65, and wavelet coefficient number of plies S is increased by 1, go to step (1) again, again operation iterative calculation
Process, PSNR > 36 of this obtained reconstructed image then stop operation, export reconstructed image.
On the basis of above-mentioned algorithm executes step, the prior-constrained MMV of structuring is introduced to of the present invention here
Restructing algorithm is illustrated compared to the reconstruction property gain that existing MMV restructing algorithm can be provided.We obtain reconstruct
The Relative reconstruction error and Y-PSNR of wavelet coefficient investigate existing 3 kinds of MMV restructing algorithms as evaluation criterion respectively
The standard that OMP-MMV, M-FOCUSS, MSBL and cM-SSBL algorithm proposed by the present invention are 1024 × 1024 to resolution ratio is surveyed
Attempt the effect being reconstructed as " cameraman ".
Fig. 4 illustrates Relative reconstruction error with the situation of change of normalization sample rate M/N, it can be seen from the figure that returning
One changes sample rate in the variation range from 0.25 to 0.55, cM-SSBL algorithm and OMP-MMV, M-FOCUSS, and MSBL algorithm is compared
It is obviously improved in terms of error performance.By taking experimental result when normalizing sample rate and being 0.4 as an example, cM-SSBL algorithm
Relative reconstruction error is 0.03 or so, and the Relative reconstruction error of OMP-MMV, M-FOCUSS, MSBL algorithm is respectively 0.075,
0.06,0.05 or so, it is 0.4 in normalization sample rate in other words, test image is that resolution ratio is 1024 × 1024
In the case where " cameraman ", cM-SSBL algorithm is compared with OMP-MMV algorithm, precision improvement 60% or so, with M-FOCUSS
Algorithm is compared, precision improvement 50% or so, and compared with MSBL algorithm, precision improvement 40% or so.
In the case that Fig. 5 is illustrated using Y-PSNR as evaluation criterion, OMP-MMV, M-FOCUSS, MSBL, and
The effect that " cameraman " test image that resolution ratio is 1024 × 1024 is reconstructed in cM-SSBL algorithm proposed by the present invention
Fruit, it can be seen from the figure that in normalization sample rate in the variation range from 0.25 to 0.55, cM-SSBL algorithm and OMP-
MMV, M-FOCUSS, MSBL algorithm are equally obviously improved compared to the aspect of performance in Y-PSNR.
In order to compare OMP-MMV, M-FOCUSS, MSBL and cM-SSBL algorithm proposed by the present invention from visual effect
Performance, we use resolution ratio to carry out for 1024 × 1024 standard testing image " cameraman " to above-mentioned three kinds of methods
Image reconstruction test, normalization sample rate are 0.4.From fig. 6, it can be seen that the quality reconstruction of cM-SSBL algorithm will be substantially better than
OMP-MMV, M-FOCUSS and MSBL algorithm.Meanwhile evaluated from objective PSNR angle, cM-SSBL method compared to
OMP-MMV, M-FOCUSS, MSBL algorithm also provide 5 gains for arriving 8dB.This explanation, introducing structuring proposed by the present invention are first
It is effective for testing the idea of constraint.
In order to examine the robustness of cM-SSBL algorithm, we have chosen different test images to OMP-MMV, M-
FOCUSS, MSBL and cM-SSBL algorithm are tested, and are counted Relative reconstruction error and Y-PSNR respectively and are carried out pair
According to experimental result is as shown in table 1.Table 1 be using different test images (resolution ratio of all test images is 1024 ×
1024, normalization sample rate image reconstructing method (cM-SSBL) proposed by the present invention and existing is based on MMV mould when being 0.4)
Property of the sparse restructing algorithm (OMP-MMV, M-FOCUSS and MSBL) of type when handling the sparse reconstruction of large-scale image
It can comparison.
Algorithms of different reconstruction property comparison of the table 1 for 1024 × 1024 dimension standard testing images
Table 1 illustrates the test result of 4 kinds of standard testing images, and 4 kinds of images are " Lena " respectively, " Boats ",
" Barbara ", and " Goldhill ", all test image resolution ratio are 1024 × 1024, and normalization sample rate is 0.4.From
Table 1 can visually see, cM-SSBL algorithm under different test images with OMP-MMV, M-FOCUSS and MSBL algorithm
It compares, can be obviously improved reconstruction property.
From figure 7 it can be seen that the number of iterations is more, the Y-PSNR of reconstruct is higher, further considers timeliness
Can, we preferably go out T in embodimentsmax=65 as the iteration total degree in emulation example.
A specific embodiment of the invention is described in conjunction with attached drawing above, but these explanations cannot be understood to limit
The scope of the present invention, protection scope of the present invention are limited by appended claims, any in the claims in the present invention base
Change on plinth is all protection scope of the present invention.
Claims (1)
1. improving sparse image reconstruction accuracy reduction complexity method with structuring is prior-constrained, which is characterized in that be a kind of right
Resolution ratio is the remote sensing images X of N × N=1024 × 10240After being redescribed using more observation vector models, X is recycledo
There is sparse coefficient on wavelet transformed domain the probability dependency of quaternary tree shape structure and more observation vector matrixes to arrange
The structuring of spatial coherence these two aspects between vector it is prior-constrained come reconstruct resolution ratio be N × N=1024 × 1024
Remote sensing images X*Method, successively realize according to the following steps in a computer:
Step (1) inputs following parameter to computer:
The remote sensing images matrix X of N × N=1024 × 10240, abbreviation matrix X0,
It is generated with the tool box MATLAB and is based on the matrix XoQuaternary tree shape structure 2-d discrete wavelet transformation of coefficient matrix
Θ, resolution ratio are N × N=1024 × 1024, and the number of plies S of two-dimensional discrete wavelet conversion is based on the matrix X0Resolution ratio
For N × N=1024 × 1024, the maximum value S of the number of pliesmax=8, level number s=1.2 ..., 8, the small echo of quaternary tree shape structure
The number of each layer of wavelet coefficient is (NC) in coefficients, root node is in first layer S1, appoint and take S=4,
With the tool box MATLAB generate M × N random observation matrix, M≤N, observing matrix abbreviation matrix Φ, similarly hereinafter,
Based on the matrix Θ generate M × N more observation vectors composition matrix Y, this be with the matrix Φ for observation
For the matrix X under conditions of matrix0Observation, i.e. Y=Φ Θ,
The remote sensing images X obtained after input reconstruct*Two image property indexs: the setting value of Y-PSNR PSNR, mean square error
The setting value of poor MSE,
Step (2) initializes Modulus Model parameter, comprising:
1. the regularized covariance Matrix C of the wavelet coefficient of different level numbers1, C2..., C4, when initialization is all unit matrix, letter
Claim Cs={ C1, C2..., C4,
2. controlling each CsThe parameter lambda of value, enables λ be initialized as 2,
3. the precision parameter for the Gaussian Profile that observation noise is obeyed is the α reciprocal of the variance of observation noisen, it is enabled to be initialized as
αn=10,
4. the precision parameter α for the Gaussian Profile that the wavelet coefficient of different level numbers is obeyed1, α2..., α4, it is each layer wavelet coefficient
The inverse of variance between column vector, when initialization, are equal to 1,
5. setting: control the probability that each layer wavelet coefficient takes zero:
For first layer S1: control matrix of wavelet coefficients Θ(1)In elementValue take zero probability to be initialized asAnd enable W1=0.3, correspondingly,The probability of negated zero is 1-W1, wherein symbol Θ(1)、In on be designated as the level number of each layer matrix of wavelet coefficients Θ, subscript (ij) is sequence of each element W in row, column
Number, similarly hereinafter,
For the second layer to the 4th layer of S2, S3, S4:
If matrix of wavelet coefficients Θ(s), the element of 2≤s≤4The coefficient value of its father's node under quad-tree structureThen controlValue takes zero probability2 s≤4 <, and initialize W0s=0.3, symbol W0sIn
Symbol " 0 " indicate father's nodeCoefficient value be zero, symbol " π " indicate father's node, similarly hereinafter,
If matrix of wavelet coefficients Θ(s), the element of 2≤s≤4The coefficient value of its father's node under quad-tree structureThen controlValue takes zero probability2 s≤4 <, and initialize W1s=0.7, symbol W1sIn
Symbol " 1 " indicate father's nodeCoefficient value be 1,
6. system hyper parameter a0, b0, c0, d0, e0, f0, the hyper parameter is the priori that the parameter in a kind of pair of step (2) is obeyed
The parameter obtained after probability distribution parameters, in which:
a0, b0, c0, d0In initialization all equal to 10-6, it is constant for controlling the parameter of probability-distribution function,
Symbol e0It include: (ec)1, and (ecz)s, 2≤s≤4:
(ec)1For controlling W1The probability density function obeyed, value are as follows:
(ec)1=0.9 × (NC)1,
(ecz)2, (ecz)3, (ecz)4It is respectively used to W0 when 2≤s of control≤42, W03, W04The probability density distribution letter obeyed
Number, value are as follows:
(ecz)2=0.1 × (NC)2, (ecz)3=0.1 × (NC)3, (ecz)4=0.1 × (NC)4,
(eco)2, (eco)3, (eco)4It is respectively used to W1 when 2≤s of control≤42, W13, W14The probability density distribution letter obeyed
Number, value are as follows:
(eco)2=0.5 × (NC)2, (eco)3=0.5 × (NC)3, (eco)4=0.5 × (NC)4,
Symbol f0It include: (fc)1, and (fcz)s, 2≤s≤4:
(fc)1For controlling W1The probability density function obeyed, value are as follows:
(fc)1=0.1 × (NC)1,
(fcz)2, (fcz)3, (fcz)4It is respectively used to W0 when 2≤s of control≤42, W03, W04The probability density distribution letter obeyed
Number, value are as follows:
(fcz)2=0.9 × (NC)2, (fcz)3=0.9 × (NC)3, (fcz)4=0.9 × (NC)4,
(fco)2, (fco)3, (fco)4It is respectively used to W1 when 2≤s of control≤42, W13, W14The probability density distribution letter obeyed
Number, value are as follows:
(fco)2=0.5 × (NC)2, (fco)3=0.5 × (NC)3, (fco)4=0.5 × (NC)4,
7. the parameter t, T of initialization control iterative processmax, in which:
T is the serial number of the number of iterations, TmaxFor maximum number of iterations, thus the setting value of the serial number of the number of iterations be t=1,
2 ..., Tmax}；
Step (3) successively calculate by following procedure the iterative process of matrix of wavelet coefficients Θ:
Step (3.1) successively carries out first time iteration according to the following steps, enables t=1,
Step (3.1.1) sets various system model parameters in constraint condition and step (2):
Constraint condition are as follows: the remote sensing images matrix X of N × N=1024 × 10240, the matrix of wavelet coefficients of N × N=1024 × 1024
Θ is quad-tree structure, Smax=8, appoint and takes S=4,
When initial, the matrix of wavelet coefficients Θ of each layer is full null matrix,
The initial value according to system model parameter and wavelet coefficient is successively calculated as follows in step (3.1.2), calculates intermediate variable
Wherein:
||·||2Indicate the 2- norm of amount of orientation or matrix,
The element Θ arranged for the i-th row of Θ matrix, jthijThe precision parameter for the probability distribution obeyed,
The element Θ arranged for the i-th row of Θ matrix, jthijThe mean of a probability distribution parameter obeyed,
The element Θ arranged for the i-th row of Θ matrix, jthijThe adjustment parameter of zero probability is taken,
αsAnd αnSystem model parameter value when expression is initial,System model parameter when expression is initial
The intermediate variable being calculated,
Φ-iIndicate to remove the matrix of other all column compositions except its i-th column in Φ,
Φ.iIndicate the vector of all elements composition of the i-th column of Φ matrix,Indicate vector Φ.iTransposition,
Θ(-i)jIt indicates in vector theta.jThe vector of the middle other all elements composition removed except its i-th of element,
Step (3.1.3), according to intermediate variable obtained in step (3.1.2)Wavelet coefficient Θ is calculated to be obeyed
Posterior probability distribution p (Θij):
Wherein: δ0The uni-impulse function at origin is represented,
Posterior probability distribution p (the Θ obeyed according to wavelet coefficient Θij), sampling obtains the i-th row, the jth column element of matrix Θ
ΘijValue to get arrive matrix Θ sampled value,
Step (3.1.4), according to the sampled value of matrix Θ obtained in the initial value of system model parameter and step (3.1.3),
Calculate the posterior probability Density Distribution and Matrix C that system model parameter is obeyed in first time iterationsFirst time iteration knot
Fruit:
p(W1)=Beta ((ec)1+||Θ(1)||0, (fc)1+N1-||Θ(1)||0),
Wherein:
The value set of footmark s is { 1 ..., 4 }, it is noted that only in the S of input >=2, and the number of plies s currently considered is greater than 1
When, just there is W0sAnd W1s,
The Frobenius norm of representing matrix,
The 0- norm of representing matrix (or vector),
Θ(s)It is that s layers of wavelet conversion coefficient combines the matrix to be formed according to natural order,
RsIt is ΘsLine number,Indicate ΘsR row,
Θs,kIndicate s layers of k-th of wavelet coefficient, Θπ(s,k)It is Θs,kFather's node coefficient value,
P (x)=Gamma (a, b) indicates that stochastic variable x obeys the gamma that parameter is (a, b) and is distributed, and expanded form isWherein Γ () indicates gamma function, by taking a is independent variable as an example, the expanded form of gamma function
For
P (x)=Beta (a, b) indicates that stochastic variable x obeys the beta that parameter is (a, b) and is distributed, and expanded form is
Π () indicates indicative function, i.e., functional value takes 1 when the expression formula in bracket is true, otherwise functional value takes 0,
According to above system model parameter αn、αs、W1、W0sAnd W1sThe posterior probability Density Distribution obeyed joins system model
Number is sampled, and updates corresponding system model parameter value, for iterative process next time,
Step (3.1.5), the sampled value of matrix Θ obtained in storing step (3.1.3), first time iterative process formally terminate,
Step (3.2) carries out second to T according to the description of step (3.1.1)~step (3.1.5)maxThe meter of secondary iteration
It calculates,
Step (3.3), to the serial number of all the number of iterations, t={ 1 ..., Tmax, calculate the institute stored in each iteration
There is the arithmetic mean of instantaneous value of matrix Θ, obtained average value is denoted as Θ*, to average value Θ*The 2-d wavelet inverse transformation for carrying out S layers, obtains
To X*,
Step (3.4), counts reconstructed image X*Relative to original image X0Error extension, PSNR and MSE, if the two
Index meets the performance indicator requirement of input, then algorithm stops, the wavelet coefficient number of plies S otherwise inputted in set-up procedure (1), and
S:=S+1 is enabled, go to step (2) again, again operation initialization and iterative process, until reconstructed image is relative to original
The error extension of beginning image meets performance indicator requirement.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710969340.9A CN107784278B (en) | 2017-10-18 | 2017-10-18 | Sparse image reconstruction accuracy reduction complexity method is improved with structuring is prior-constrained |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710969340.9A CN107784278B (en) | 2017-10-18 | 2017-10-18 | Sparse image reconstruction accuracy reduction complexity method is improved with structuring is prior-constrained |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107784278A CN107784278A (en) | 2018-03-09 |
| CN107784278B true CN107784278B (en) | 2019-02-05 |
Family
ID=61433781
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710969340.9A Active CN107784278B (en) | 2017-10-18 | 2017-10-18 | Sparse image reconstruction accuracy reduction complexity method is improved with structuring is prior-constrained |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107784278B (en) |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102722866A (en) * | 2012-05-22 | 2012-10-10 | 西安电子科技大学 | Compressive sensing method based on principal component analysis |
| CN103093445A (en) * | 2013-01-17 | 2013-05-08 | 西安电子科技大学 | Unified feature space image super-resolution reconstruction method based on joint sparse constraint |
| CN104657944A (en) * | 2014-12-31 | 2015-05-27 | 中国科学院遥感与数字地球研究所 | Compressed sensing remote sensing image reconstruction method based on reference image texture constraints |
| CN106204666A (en) * | 2015-06-12 | 2016-12-07 | 江苏大学 | A kind of compression sensed image reconstructing method |
| WO2017125926A2 (en) * | 2016-01-18 | 2017-07-27 | Dentlytec G.P.L. Ltd | Intraoral scanner |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102148987B (en) * | 2011-04-11 | 2012-12-12 | 西安电子科技大学 | Compressed sensing image reconstructing method based on prior model and 10 norms |
-
2017
- 2017-10-18 CN CN201710969340.9A patent/CN107784278B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102722866A (en) * | 2012-05-22 | 2012-10-10 | 西安电子科技大学 | Compressive sensing method based on principal component analysis |
| CN103093445A (en) * | 2013-01-17 | 2013-05-08 | 西安电子科技大学 | Unified feature space image super-resolution reconstruction method based on joint sparse constraint |
| CN104657944A (en) * | 2014-12-31 | 2015-05-27 | 中国科学院遥感与数字地球研究所 | Compressed sensing remote sensing image reconstruction method based on reference image texture constraints |
| CN106204666A (en) * | 2015-06-12 | 2016-12-07 | 江苏大学 | A kind of compression sensed image reconstructing method |
| WO2017125926A2 (en) * | 2016-01-18 | 2017-07-27 | Dentlytec G.P.L. Ltd | Intraoral scanner |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107784278A (en) | 2018-03-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Metzler et al. | Unsupervised learning with Stein's unbiased risk estimator | |
| CN109767007B (en) | Minimum mean square error detection method based on quantum computation | |
| CN109116293B (en) | Direction-of-arrival estimation method based on lattice-separated sparse Bayes | |
| CN107817465A (en) | The DOA estimation method based on mesh free compressed sensing under super-Gaussian noise background | |
| CN109298383B (en) | Mutual-prime array direction-of-arrival estimation method based on variational Bayes inference | |
| CN110109050A (en) | The DOA estimation method of unknown mutual coupling under nested array based on sparse Bayesian | |
| CN105761223A (en) | Iterative noise reduction method based on image low-rank performance | |
| Urschel et al. | A cascadic multigrid algorithm for computing the Fiedler vector of graph Laplacians | |
| Cosme et al. | Implementation of a reduced rank square-root smoother for high resolution ocean data assimilation | |
| CN114786018A (en) | Image reconstruction method based on greedy random sparse Kaczmarz | |
| CN110174658B (en) | Direction-of-arrival estimation method based on rank-dimension reduction model and matrix completion | |
| CN115453528A (en) | Method and device for realizing segmented observation ISAR high-resolution imaging based on rapid SBL algorithm | |
| Jensen et al. | Simulation of conditioned semimartingales on riemannian manifolds | |
| CN107784278B (en) | Sparse image reconstruction accuracy reduction complexity method is improved with structuring is prior-constrained | |
| CN112087235A (en) | Sparsity self-adaptive DOA estimation method and system based on pseudo-inverse perception dictionary | |
| CN110174657B (en) | Direction-of-arrival estimation method based on rank-one dimension reduction model and block matrix recovery | |
| CN114972041B (en) | Polarization radar image super-resolution reconstruction method and device based on residual error network | |
| Wang et al. | Image reconstruction from patch compressive sensing measurements | |
| CN108896967B (en) | Method and device for detecting distance extended target based on clutter covariance matrix estimation | |
| Gilman et al. | Online tensor completion and free submodule tracking with the T-SVD | |
| CN114844544A (en) | Low-tube rank tensor decomposition-based co-prime array beamforming method, system and medium | |
| CN113838104B (en) | Registration method based on multispectral and multimodal image consistency enhancement network | |
| Gunn et al. | Regularized training of intermediate layers for generative models for inverse problems | |
| Pu et al. | ORTP: A video SAR imaging algorithm based on low-tubal-rank tensor recovery | |
| He et al. | TDiffDe: A Truncated Diffusion Model for Remote Sensing Hyperspectral Image Denoising |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |