fast acquisition is an important issue in magnetic resonance imaging ( mri ) for avoiding physiological effects and reducing scanning time on patients .
unfortunately , accelerating acquisition by reducing k - space samples leads to noise amplification , blurred object edges , and aliasing artifacts in mr reconstructions .
then improving reconstruction accuracy from highly undersampled k - space data becomes a complementary tool to alleviate the above side effects of reducing acquisition . to cope with the loss of image quality
one process of introducing prior information in the reconstruction is known as  regularization  [ 211 ] .
tikhonov regularization , a commonly used method that pursuits reconstructions by a l 
 2-norm minimization , leads to a closed - form solution that can be numerically implemented in an efficient way [ 68 , 11 ] .
moreover , with the advent of compressed sensing ( cs ) theory , sparsity - promoting regularization has gained popularity in mri ( e.g. , the l 
 1-based regularization ) [ 13 , 5 , 10 , 12 ] . the cs theory states that sparse , or more generally , compressible signals , incoherently acquired in an appropriate sense , can be recovered from a reduced set of measurements that are largely below the nyquist sampling rate .
exact reconstruction can be achieved by nonlinear algorithms , using l 
 1 minimization or orthogonal matching pursuit ( omp ) [ 13 , 5 , 10 , 12 , 13 ] .
< p < 1 ) minimization will guarantee a better recovery by directly attacking the l 
 0 minimization problem [ 1416 ] .
1 ) norm optimization as a relaxation of the l 
 0-norm for mri reconstruction [ 17 , 18 ] and showed that it often outperforms the l 
 1-minimization method .
inspired by viewing p as a scale parameter , trzasko and manduca   developed a homotopic l 
 0-minimization strategy to reconstruct mr image and generalized the nonconvex l 
 p ( 0
cands et al . presented the iteratively reweighted l 
 1 minimization ( irl1 ) algorithm , which amounts to linearly minimize the log penalty of the total variation .
wong et al .   further incorporated a seminonlocal priority into the homotopic l 
 0-minimization for breast mri . besides the prior - inducing penalty function , another important issue in cs - mri reconstruction , is the choice of a sparsifying transform , since one of the important requirements in cs - mri is that the image to be reconstructed has a sparse representation in a given transform domain . at the start , sparse representation was usually realized by total variation ( tv ) regularization and/or wavelet transform [ 1 , 12 , 13 , 15 , 20 , 21 ] .
as it is known , tv prior assumes the target images consist of piecewise constant areas which may not be valid in many practical mri applications .
other analytically designed bases , such as wavelets and shearlets , also involve their intrinsic deficiencies .
knoll et al . introduced the total generalized variation ( tgv ) for mr image reconstruction problems .
liang et al .   applied the nonlocal total variation ( nltv ) regularization to improve the signal - to - noise ratio ( snr ) in parallel mr imaging , by replacing the gradient functional in conventional tv using a weighted nonlocal gradient function to reduce the blocky effect of tv regularization .
since fixed bases might not be universally optimal for all images , some methods utilizing sparsity prior under adaptive transform / dictionary were developed recently [ 2429 ] .
the sparsity in this framework is enforced on overlapping image patches to emphasize local structures .
additionally , the global dictionary consisting of local atoms ( prototype ) is adapted to the particular image instance , which thereby provides better sparsity .
our work also adopts the similar concept as a penalty term for mri reconstruction . after reviewing two important components in cs
, an interesting question is whether a better sparsity - inducing function ( e.g. , l 
 p ( 0
1 ) norm ) combined with a better transform ( e.g. , adaptive dictionary ) will lead to better reconstruction from the same set of measurements . to the best of our knowledge , there is only one paper by shi et al . preliminarily studying the l 
 0-approximation sparsity under the learned dictionary , where the l 
 0 penalty is relaxed by the nonconvex minimax concave ( mc ) penalty .
however , their work focused on image processing and the homotopic strategy they used is computationally demanding . in this paper
, we propose a novel method to solve the nonconvex minimization problem in cs - mri by incorporating the reweighting scheme with our recently developed two - level bregman method with dictionary updating ( tbmdu ) .
an ultimate advantage of our proposed method is that , by incorporating the iteratively reweighted scheme into the two - level bregman method / augmented lagrangian ( al ) with dictionary updating , the l 
 p - based updating step is conducted in each inner iteration of the al scheme , and both the minimization of the sparsity of image patches and the data - fidelity constraint can be solved in a unified formalism . the modified alternating direction method ( adm ) efficiently solves the model of pursuing the approximated l 
 p ( 0 < p < 1)-norm penalty . specifically , by applying iteratively reweighted l 
 2 or l 
 1 minimization scheme at each al inner iteration
, the proposed algorithm utilizes adaptive weights to encourage the large coefficients to be nonzero and the small ones to be zero under learned dictionary , hence resulting in more compact and discriminative representation for reconstructed image patches .
numerical experiments show that the performance of the derived method is superior to other existing methods under a variety of sampling trajectories and k - space acceleration factors .
particularly , it achieves better reconstruction results than those by using tbmdu without adding any computational complexity .
we start with a brief review on tbmdu and iteratively reweighting scheme for l 
 p - minimization in section 2 .
several numerical simulation results are illustrated in section 4 to show the superiority of the proposed method .
plenty of papers , with the theme of compressed sensing ( cs ) [ 2 , 3 ] , demonstrated that mr images can be reconstructed from very few linear measurements [ 13 , 5 , 10 , 12 ] .
the well - known example is the success of tv regularization which implies that mr images can be approximated by those having sparse gradients .
generally , if a sparsifying transform  can be readily defined , ideally one could choose the estimation with the sparsest representation in  that still matches the limited measurements .
that is , by minimizing transform - domain sparsity - promoting energy j(u ) = ||u||0 that is subject to the data consistency f = f 
 p 
 u + n , where n is noise , it yields 
 ( 1)minu  j(u)s.t . 
||fpuf||2 , 
 where u  c 
  is a vector representing the n - pixel 2d complex image .
f 
 p  c 
  denotes the partially sampled fourier encoding matrix and f  c 
  represents the acquired data in k - space .
=2 is the standard deviation of the zero - mean complex gaussian noise , and  is the standard deviation of both real and imaginary parts of the noise .
practically , the l 
 0 quasi - norm ||||0 in the regularization term of ( 1 ) is usually relaxed . in this paper , we propose some efficient algorithms for nonconvex relaxation such as l 
 p quasi - norm ||||p , 0 < p < 1 [ 17 , 18 ] . the choice of the sparsifying transform  is an important consideration in minimizing functional ( 1 ) . besides the fixed transforms such as finite - difference and wavelet transform , modeling image patches ( signals ) by sparse and redundant representations have been drawing considerable attention in recent years . considering the image patches set ru = [ r 
 1 
 u , r 
 2 
 u ,  ,
r 
 l 
 u ] consisting of l signals , r 
 l 
 u  c 
  denotes a vectored form of the mm patch extracted from the image u of size nn .
the sparseland model for image patches suggests that each image patch , r 
 l 
 u , could be sparsely represented over a learned dictionary d 
  ; that is , 
 ( 2)l=argminl(||l||0+2||dlrlu||22 ) ,  l=1,2,  ,l , 
 where  
 l stands for the sparse coefficient of the lth patch .
the combination of sparse and redundant representation modeling of signals , together with a learned dictionary using signal examples , has shown its promise in a series of image processing applications [ 2429 ] .
for example , in [ 25 , 26 ] , the authors used reference mri image slices to train the dictionary and achieved limited improvement . considering that a dictionary learnt from a reference image would not be able to effectively sparsify new features in the current scan ; ravishankar and bresler   suggested learning dictionary from the target image itself and developed a two - step alternative method to remove aliasing artifacts and noise in one step and subsequently fills in the k - space data in the other step . in the following , to make the paper self - contained , we , respectively , review the tbmdu method for mri reconstruction and reweighted scheme for solving nonconvex relaxation of the l 
 p - quasi - norm minimization .
recently , we proposed a series of dictionary learning methods [ 3032 ] based on augmented lagrangian / bregman iterative scheme . particularly , by using a two - level bregman iterative method and the relaxed version of sparseland model ( i.e. , j(u ) = l(|| 
 l||1 + ( /2)||d 
 lr 
 l 
 u||2 
 ) ) as the regularization term for mri reconstruction , the tbmdu method solves the objective function equation ( 1 ) by iterating the following : 
 ( 3)uk+1=argminu{mind,l(||l||1+2||dlrlu||22 )  + 2||fpufk||22},fk+1=fk+ffpuk+1 , 
 where d = [ d 
 1 , d 
 2 ,  , d 
 j ]  c 
  and  = [  
 1 ,  
 2 ,  ,  
 l ]  c 
 .
 stands for the sparse level of the image patches in the  optimal  dictionary . for many natural or medical images
, the value of  can be determined empirically with robust performance in our work .
j = k  m , in which k measures the degree of the overcompleteness of the dictionary .
the proposed tbmdu method consists of a two - level bregman iterative procedure , where ( 3 ) is the outer - level of the method . on the other hand ,
the inner - level bregman iterative is employed to solve the subproblem of sparse representation in ( 3 ) , that is , adding auxiliary variables z 
 l to convert the unconstrained problem ( 2 ) into a constrained formulation : 
 ( 4)mind,l  ||l||1+2||zl||22s.t . 
then split bregman method / augmented lagrangian is used to solve problem ( 4 ) as follows : 
 ( 5)(di+1,li+1,zli+1 )  = argmind,l , zll||l||1+2||zl||22+2||zl+dlrluyli||22 , 
 ( 6)yi+1=yi+(ruzi+1di+1i+1 ) . 
 
we minimize the corresponding al function ( 5 ) alternatively with respect to one variable at a time .
this scheme decouples the minimization process and simplifies the optimization task . at the dictionary updating step , by updating ( 5 ) with respect to d along the gradient descent direction and constraining the norm of each atom to be unit
, we can get the following update rule : 
 ( 7)di+1=di+yi+1(i+1)t , dji+1=dji+1||dji+1||2 ,  j=1,2,  ,j . 
  at the sparse coding stage of updating  
 l ,
the optimization problem for each  
 l is derived : 
 ( 8)li , m+1=argminl||dilblyli||22 + 2+||l||1 . 
  by applying the iterative shrinkage / thresholding algorithm ( ista ) [ 33 , 34 ] ( i.e. , performing a gradient descent of the first functional term in ( 8) and then solving a proximal operator problem ) and using the immediate variable y 
 l 
  to represent the updating scheme compactly ( i.e. , the formulation y 
 l 
  = ( / + )(d 
  
 l 
  + b 
 l + y 
 l 
 / ) derived from ( 6 ) )
, it attains the solution of  
 l as follows : 
 ( 9)li , m+1=shrink(li , m+(+)(di)tylm+1,(+) ) , 
 where  > eig((d 
 ) 
 d 
 ) and x = shrink(g ,  ) = ( g  max(|g|   , 0))/(max(|g|   , 0 ) +  ) is the solution of x = arg  minx||x||1 + ( 1/2)||xg||2 
 .
line 10 is the frequency interpolation step for updating the image u. al scheme and adm applied to the constrained problem result in an efficient two - step iterative mechanism , which alternatively updates the solution u and image - patch related coefficients (  , y , and d ) . for a more detailed description of the derivations , the interested readers can refer to . in ( 1 )
for example , some nonconvex relaxations such as l 
 p quasi - norm ||||p , 0 < p < 1 [ 17 , 18 ] , log penalty , and smoothly clipped absolute deviation ( scad ) penalty   were usually used .
compared to l 
 1 norm , l 
 p function with 0 < p < 1 is a smooth relaxation of l 
 0 and the geometry of l 
 p ball gives better approximation on sparsity .
recent numerical results showed that adapting nonconvex optimization technique can reduce the required number of measurements for reconstruction [ 1416 ] .
most current approaches for solving nonconvex l 
 p - quasi - norm optimization can be classified into two categories : reweighted l 
 2 and reweighted l 
 1[3641 ] .
the common idea of these methods is to utilize an optimization transfer technique to iteratively surrogate the nonconvex potential function by a convex function [ 37 , 39 , 41 ] .
usually , the convex function is chosen as a local quadratic approximation or linear local approximation .
a technique drawing much interest is the iteratively reweighted least squares ( irls ) method which is widely used in compressed sensing and image processing [ 36 , 40 ] .
another technique closely related to irls is the iteratively reweighted norm ( irn ) approach introduced by rodrguez and wohlberg for generalized total variation functional .
the irn method pursues to minimize the l 
 p norm f( ) = ( 1/p)||||p 
  for p  2 by using a weighted l 
 2 norm in an iterative manner . at iteration
k + 1 , the solution  
  is the minimizer of 
 ( 10)g()=f~(,k)=12||wk||2 ,  wk = diag(|k|p2 ) . 

the sequence of solution {  
 } converges to the minimizer of f( ) as k  + . it has been proven that irn is one kind of majorization - minimization ( mm ) method , which involves good property of g( 
 )  g( 
 ) . in order to extend the iterative reweighting approach irls to the general nonconvex function , mourad and reilly developed the quadratic local approximation of the nonconvex function as follows : 
 ( 11)f(j , l)=(|j , l|)=(|j , lk|)+(|j , lk|)(|j , l||j , lk|)+0(|j , l||j , lk|)2 , 
 where  
 0  0 is a suitably chosen parameter and  is the first derivative of . j stands for the jth index number of vector  
 l. a similar work conducted by elad and aharon is the iteratively reweighted l 
 1 minimization ( irl1 ) algorithm , which iteratively minimizes the linearization of a quasi - logarithmic potential function .
zou and li   presented the linear local approximation of the general nonconvex function and the upper bound of 
 ( 12)f(j , l)=(|j , l|)=(|j , lk|)+(|j , lk|)(|j , l||j , lk| ) . 

equation ( 12 ) aims to approximate the penalty function with its first - order taylor expansion at the current iteration .
its majorization properties can be easily obtained by exploiting the concavity of the function  , which always lies below its tangent .
the numerical tests in   by chartrand and yin demonstrated that irls is comparable with irl1 . in , the authors presented weighted nonlinear filter for compressed sensing ( wnfcs ) for general nonconvex minimization by integrating the irl1 algorithm shown in ( 12 ) and their previously developed nfcs framework .
as discussed in previous section , dictionary updating and sparse coding to ( 5 ) are performed sequentially .
therefore , an interesting question is that will the better sparsity inducing function in the coefficient matrix lead to better image reconstruction under the learned dictionary ? in the following , we investigate two variations at the sparse coding step . by employing the l 
 p - seminorms ( 0 < p  1 ) , it is feasible to expect that utilizing approximations that are closer to the l 
 0-quasi - norm than l 
 1 will correspondingly reduce the required number of measurements for accurate reconstruction . by replacing l 
 1 norm with nonconvex
l 
 p norm , ( 3 ) can be rewritten as 
 ( 13)uk+1=argminu{mind,l(1p||l||p+2||dlrlu||22 )  + 2||fpufk||22}.fk+1=fk+ffpuk+1 . 
  after some manipulations similar to those in section 2.1 , updating coefficients at the sparse coding stage are deduced as follows : 
 ( 14)li , m+1=argminl||dilblyli||22 + 2+1p||l||p . 
 
when other variables are fixed , numerical strategies are proposed in the following to solve the corresponding nonconvex minimization problems with regard to  , by iteratively minimizing a convex majorization of the nonconvex objective function . by inserting the local linear approximation in the sparse coding step ,
it is possible to transform the nonconvex unconstrained minimization problem into a convex unconstrained problem , which can be solved iteratively . in the penalty functional of ( 14 ) , like the well - known method ista , a proximal operator is employed to approximate the first term , and the reweighted l 
 1-approximation is employed to the second term : 
 ( 15)li , m+1=argminl||dilblyli||22ista+2+1pj|j , l|prel1 , 
 where j stands for the jth index number of vector  
 l. following the definition in ( 12 ) , where (| 
 j , l| ) = ( 1/p)| 
 j , l| and (| 
 j , l| ) = 1/(| 
 j , l| +  ) , it attains the solution of  
 l : 
 ( 16)li , m+1=argminl  ||l[li , m+(+)(di)tylm+1]||22  + 2+||wli , ml||1=shrink(li , m+(+)(di)tylm+1,(+)wli , m ) , 
 where w 
 j , l 
  = 1/(| 
 j , l 
 | +  ) , j = 1,2 ,  , j.  > 0 is a small parameter to prevent division by zeros .
compared to the previous reweighted l 
 1 strategy , the only difference is that the reweighted l 
 2-approximation in ( 10 ) is employed to the second term : 
 ( 17)li , m+1=argminl||dilblyli||22ista+2+1p||l||pprel2 . 
 
then it attains the solution of  
 l : 
 ( 18)li , m+1=argminl  ||l[li , m+(+)(di)tylm+1]||22++||wli , ml||22 . 
 
the minimizer of this least - square problem is given by 
 ( 19)li , m+1=li , m+(+)(di)tylm+1/+(+)wli , m/ , 
 where w 
 j , l 
  = 1/(| 
 j , l 
 | +  ) , j = 1,2 ,  , j. similar to that in ( 14 ) ,  > 0 is a small parameter to prevent numerical instabilities .
now , we summarize our proposed method for mri reconstruction here , which we call weighted tbmdu . the detailed description of the proposed method is listed in algorithm 2 .
similar to tbmdu , the proposed wtbmdu method alternatively updates the target solution u , image patch related coefficients (  , y , and d ) , and auxiliary variables .
the difference between the plain tbmdu and the weighted tbmdu mainly lies on the weights used in updating the sparse coefficients . in wtbmdu ,
the weights are obtained by evaluating the function of variables at the solution of the previous step .
specifically , whether the w 
 j , l 
  = 1/(| 
 j , l 
 | +  ) in ( 16 ) or w 
 j , l 
  = 1/(| 
 j , l 
 | +  ) in ( 19 ) , the weighting scheme w 
 j , l 
  decreases as the absolute value of  
 j , l 
  increases , indicating that it penalizes more on the coefficients with small magnitude value .
therefore , this operation strongly encourages the large coefficients to be nonzero and the small ones to be zero . in line 8 for reweighted l 
 1 and line 10 for reweighted l 
 2 ,
the formulations denote operating every element in the matrix component - wise . in the following content of this paper , we denote the reweighted l 
 1 and reweighted l 
 2 in wtbmdu as wtbmdu - l1 and wtbmdu - l2 , respectively .
the strategy we used in wtbmdu is similar to that used in wnfcs , that is , combining the penalized proximal splitting strategy and reweighting strategy .
the difference is that wnfcs focus on updating u in the simple cs domain and the proximal splitting strategy is employed to the fourier - related data - fidelity term , while our wtbmdu devotes to updating the coefficients in the adaptive dictionary and the proximal splitting strategy is employed to the dictionary - related sparse representation error term .
the proposed method involves four parameters :  ,  ,  , and . the setting of these parameters is similar to that in tbmdu .
firstly , both parameters  and  are the bregman ( or augmented lagrangian ) positive parameters associated with the small image patches and the whole image itself , respectively .
one is for the overlapping image patches and the other is for the image solution itself .
it has been mathematically proven that the choice of the bregman parameter has little effect on the final reconstruction quality as long as it is sufficiently small [ 31 , 44 , 45 ] . in our work ,
the smaller the value of the bregman parameter is , the more iterations the bregman method need to reach the stopping condition .
moreover , since  and  are with different orders of magnitude , we set   n/m in the al formalism for the balance between various 
secondly ,  stands for the sparse level of the image patches and can be determined empirically .
finally , the step size  in the dictionary updating stage can be set to be a small positive number , for example , 0.01 . in summary
, there is only the parameter  that should be carefully chosen in order to enable the efficiency of the algorithm . in our method , by taking advantage of the typical structure of the problem in this paper , we propose a similar rule to that used in .
specifically , we set  and  so as to achieve condition numbers (1 + w 
 l 
 /t 
 0 ) and (ff 
 p 
 f 
 p 
 f 
  + ( /)fl 
 r 
 l 
 r 
 l 
 f 
 ) that result in fast convergence of the algorithm . since t 
 0 = /( +  ) and    , t 
 0  . consequently (1 + w 
 l 
 /t 
 0 ) is a decreasing function of . choosing  such that (1 + w 

l 
 /t 
 0 )  1 would require a large  and accordingly , the influence of the weight w 
 l 
  diminishes and the solution of  
 l 
  trends to the least - square solution ( q = 2 ) . in our experiments , we observed that this phenomenon would result in an  immature  dictionary that is not learned enough . on the other hand , taking   0 would increase (1 + w 

l 
 /t 
 0 ) which makes ( 1 + w 
 l 
 /t 
 0 ) numerically unstable .
the same trend also applies to (ff 
 p 
 f 
 p 
 f 
  + ( /)fl 
 r 
 l 
 r 
 l 
 f 
 ) as a function of /. we found that empirically choosing  such that e((w 
 l 
 ) )  l 
 (w 
 l 
 ) /l
additionally , we can estimate  in a more practical way . specifically , since the number of data samples is huge , we can set e(w 
 l 
 )  t 
 0 , where e(w 
 l 
 ) [0.2,1 ] is a robust estimate . together with the fact that   100 , it leads to   [ 1/500 , 1/100 ] . as for the parameter  , since (ff 
 p 
 f 
 p 
 f 
  + ( /)fl 
 r 
 l 
 r 
 l 
 f 
 ) = (  + )/ , where  = 64 when the patch overlapping is r = 1 , it concludes that setting   n/m such that (ff 
 p 
 f 
 p 
 f 
  + ( /)fl 

r 
 l 
 r 
 l 
 f 
 ) = 64 + 1    is a reasonable choice , under the assumption that the k - space data is significantly undersampled ( e.g. ,   [ 5% , 30% ] ) . as for the convergence , because of the unconvexity and nonlinearity of the problem in the case of updating dictionary , the global solution may not to be found easily like in tbmdu .
nevertheless , our dictionary is updated by a gradient descent of the al scheme which leads to a monotonic decrease in the cost function . at the sparse coding stage , since wtbmdu merges the reweighted and penalization strategies , it still maintains the descent property of the latter and local minimum must be attained as demonstrated in .
therefore , both the value of the objective function and the norm of the reconstruction difference between successive iterations can be chosen as the stopping criterion .
in this section , we evaluate the performance of the proposed method using a variety of sampling schemes , with different undersampling factors .
sampling schemes used in our experiments include 2d random sampling , cartesian approximation of multishot variable - density spiral sampling , cartesian sampling with random phase encodings ( 1d random ) [ 1 , 15 ] , and pseudo radial sampling [ 15 , 29 ] .
reconstruction results on simulated mri data , a complex phantom , and real mri data were presented .
the mr images tested in the synthetic experiments are from in vivo mr scans of size 512  512 ( many of which are courtesy ( 2009 , american radiology services [ online ] .
available : http://www3.americanradiology.com/pls/web1/wwimggal.vmg/ ) and used in ) , and the real mri data examples reported here are of size 256  256 ( except in figures 6 and 7 where a phantom of size 512  512 was used ) . according to many prior work on cs - mri [ 1 , 13 , 29 ] , the cs data acquisition was simulated by subsampling the 2d discrete fourier transform of the mr images ( except in the second subsection where real acquired data was used ) .
our proposed method wtbmdu was compared with the leading dlmri ( the code is available in https://netfiles.uiuc.edu/ravisha3/www/dlmricode.zip )   and tbmdu methods , which have been shown to substantially outperform other cs - mri methods such as ldp ( matlab codes are available in http://www.stanford.edu/~mlustig/ ) , and the zero - filling reconstruction .
dlmri directly solves the l 
 0-minimization by omp while tbmdu is devoted to the l 
 1-induced sparse minimization . in
each given example , the parameters for the dlmri method were set to be default values .
in the experiments , the nominal values of various parameters were set as patch size m=6 , the over - completeness of the dictionary k = 1 ( correspondingly j = 36 ) , and the patch overlap r = 1 ; thereby the number of data samples were l = 65536 for n=256 and
=(m / r)2=36 ,  = 0.0056 ,  = n/m ,  = 0.01 ,  = 12 , and m = 3 . to avoid dividing by zero ,
the parameter  used in the weighted matrix w 
 q , l 
  decreased by 2% after each inner iteration with an initial value of 5 .
the setting of parameters m and  is very similar to that in , [ 31 , 32 ] and not discussed here due to the limit of paper space .
real - valued dictionaries were used for the simulated experiments with real - valued images , where the over - complete discrete cosine transform ( dct ) was chosen as the initial dictionary .
complex - valued dictionaries were used for real mr data , where both the real and imaginary parts were the same dct matrix .
the quality of the reconstruction was quantified using the peak signal - to - noise ratio ( psnr ( the psnr is defined as psnr = 20  log10  255/rmse , where the rmse is the root mean error estimated between the ground truth and the reconstructed image ) ) and high - frequency error norm ( hfen ) .
all algorithms were implemented in matlab 7.1 on a pc equipped with amd 2.31  ghz cpu and 3 gbyte ram .
figure 1 involves an axial t2-weighted reference image of the brain with pseudo radial sampling under 85% and 95% undersampling percentages , respectively ( i.e. , only acquiring 15% and 5% k - space data with corresponding acceleration factors of 6.67 and 20 ) . the plots of psnr and hfen values as functions of iteration number under the undersampling percentage of 85% are presented in figures 1(d ) and 1(e ) .
it can be observed that the quantitative measures of both tbmdu and wtbmdu change quickly during the first few iterations . in other words ,
these measure values only need less iterations to reach the convergence zone and hence the iterative convergence property of our method is better than that of dlmri .
the higher psnr values and lower hfen values after convergence also confirm the superiority of our method to dlmri and tbmdu .
the reconstructed results shown in figures 1(f ) , 1(g ) , 1(h ) , and 1(i ) reveal that the method wtbmdu - l1 with p = 0.5 provides a more accurate reconstruction on image contrast and sharper anatomical depiction .
compared to dlmri , the magnitude images of the reconstruction error shown in figures 1(h ) and 1(i ) indicate that our method exhibits crisper reconstruction of object edges ( the large anatomical structure in the middle region ) and preserves finer texture information ( the gray matter regions in the bottom - right of the reconstruction ) . in general , our proposed method provides better intensity fidelity to the fully sampled image .
more obvious differences in visual quality can be observed in the case of 95% undersampling as shown in figures 1(j ) and 1(k ) .
the obtained psnrs of dlmri and wtbmdu - l1 with p = 0.5 are 26.30  db and 28.13  db , respectively .
the dlmri reconstruction in figure 1(j ) based on k - svd dictionary updating and greedy pursuit of coefficients shows a large number of spurious oscillations , although it gives much improvement than zero - filling and ldp ( not shown in this paper ) .
in contrast , wtbmdu shown in figure 1(k ) results in much fewer spurious oscillations and better preservation of edges through iterative weighting of the coefficients under the data - adaptive dictionary .
figure 2 compares the results generated by dlmri and wtbmdu using four sampling trajectories roughly under the same undersampling percentage : variable density random with 87% undersampling , cartesian approximation of multishot spiral with 86% undersampling , 1d cartesian trajectory , and pseudo radial sampling with 86% undersampling ( i.e. , 7.11-fold acceleration ) .
the test image is the axial t2-weighted brain image shown in figure 1(a ) . as can be observed from the error images , wtbmdu - l1 performs better in reducing aliasing artifacts and maintaining fine details than dlmri with all verified sampling trajectories . 

table 1 lists the psnr values of the axial t2-weighted brain image at different sampling trajectories with the same undersampling percentage using dlmri , tbmdu , and wtbmdu .
generally , the improvements gained by wtbmdu over other methods are different for four kinds of trajectories although under the same undersampling rate .
the largest and smallest improvements were achieved with the 2d random and radial sampling , respectively , where roughly 5  db and 2  db were obtained .
this indicates that the efficiency of dictionary learning methods may depend on the incoherence of the data acquisition . in the family of wtbmdu algorithms ,
the optimal p value is also different at various trajectories for both wtbmdu - l1 and wtbmdu - l2 . to balance the numerical calculation and
the selection of trajectories , p = 0.5 or p = 0.7 is a good option . 

figure 3 illustrates the performance of dlmri , tbmdu , and wtbmdu at a range of acceleration factors including 2.5 , 4 , 6 , 8 , 10 , and 20 , where zero - mean complex white gaussian noise with standard deviation  = 10.2 was added to the 2d random sampled k - space . since the stopping rule for the outer loop of both tbmdu and wtbmdu is determined by ||f 
 p 
 u 

f||2 <  , the number of outer iterations k 
 max of tbmdu and wtbmdu - l2 with p = 0.7 take the values of 3 , 3 , 6 , 7 , 7 , 11 and 4 , 8 , 8 , 10 , 11 , 11 for the above six acceleration factors , respectively . for the quantitative comparison , the values of psnr and hfen as functions of acceleration factor are shown in figures 3(b ) and 3(c ) . it can be seen that wtbmdu performs best at all tested acceleration factors .
the reconstruction results and the corresponding error images under 10-fold acceleration using three methods are displayed in figures 3(d ) , 3(e ) , and 3(f ) and 3(g ) , 3(h ) , and 3(i ) , respectively .
as shown in the error images , the brighter the error image appears , the larger the deviation between the reconstruction and the reference image will be .
wtbmdu ( figure 3(i ) ) presents less pixel errors and structure loss than that of dlmri in figure 3(g ) and tbmdu in figure 3(h ) , especially in regions indicated by red arrows .
in general , it can be concluded that both wtbmdu - l2 and tbmdu offer strong preservation of details compared with that of dlmri , while wtbmdu - l2 provides a crisper result of fine structures than that of tbmdu . 

figure 4 depicts the reconstruction results of a transverse slice of a noncontrast mr angiography ( mra ) of the circle of willis ( cow ) at the same experiment setting as in figure 3 .
the mra of the circle of willis has much textural information such as the vessels on the middle region and fine - scale details on the bottom region .
as can be seen from figure 4(b ) , when the acceleration factor increased until 10-fold , the psnr gap between tbmdu and wtbmdu increases synchronously .
the psnr improvement is also reflected in figures 4(e ) and 4(f ) where much less errors appeared in the wtbmdu error image , indicating wtbmdu performs better in maintaining fine details . on the other hand , similar as observed in figure 3 , when the acceleration factor increased to as large as 20-fold , the advantage of the nonconvex optimization degraded and none of these methods can faithfully reconstruct the original image from such fewer k - space samples with additional noise . to investigate the sensitivity of various methods to different levels of complex white gaussian noise , dlmri , tbmdu , and wtbmdu
were applied to reconstruct a t2-weighted sagittal view of the lumbar spine under pseudo radial sampling at 6.09-fold acceleration .  
figure 5(c ) plots the psnr values of the recovered mr images by dlmri ( blue curves ) , tbmdu ( green curves ) , and wtbmdu ( red curves ) at a sequence of different standard deviations (  = 2 , 5 , 8 , 10 , 14.2 ) . in the case of 
= 2 , the psnr of the image obtained by dlmri is only 38.22  db , tbmdu is 39.49  db , and wtbmdu - l2 with p = 0.7 reaches 40.52  db .
the reconstructions and corresponding magnitudes of the error images with  = 8 are shown in figures 5(d ) , 5(e ) , 5(f ) , 5(g ) , 5(h ) , and 5(i ) .
it can be observed that the skeletons in the top half part of the tbmdu reconstruction appear less obscured than those in the dlmri results .
meanwhile , the reconstruction by wtbmdu is clearer and sharper than that by dlmri and tbmdu and is relatively devoid of aliasing artifacts .
this reveals that our method provides a more accurate reconstruction of image contrast and sharper anatomical depiction in noisy case . 

figure 6 shows the comparison between different methods on a physical phantom , which is often used to assess the resolution of mr reconstruction .  
figures 6(c ) , 6(d ) , 6(e ) , and 6(f ) exhibit the results of dlmri , tbmdu , wtbmdu - l1 , and wtbmdu - l2 with p = 0.7 at 80% undersampling ratio .
the corresponding psnr values are 18.66  db , 26.82  db , 28.89  db , and 29.12  db , respectively .
we can find that the wtbmdu reconstructions exhibit higher resolution than those with dlmri and tbmdu and are almost devoid of aliasing artifacts especially in zoomed - in regions . to investigate the noise sensitivity of proposed method on the complex k - space data , complex gaussian noise of  = 30
was added to the k - space with 5-fold acceleration . in this case , the psnr values of dlmri , tbmdu , and wtbmdu - l2 with p = 0.7 methods are 17.93  db , 22.94  db , and 25.42  db , respectively .
the enlargements of two region - of - interests ( rois ) are presented in figures 7(d ) and 7(e ) . as can be observed ,
the dlmri reconstruction exhibits more oscillating artifacts than that of the other two methods . besides , the spot and circle in the bottom right of the wtbmdu reconstruction appear less obscured than those in the dlmri and tbmdu results . in figures 8 and 9 ,
two real brain data sets containing more fine - detailed structures [ 44 , 46 ] were used for comparison .
the images in figures 8(a ) and 9(a ) are both the fully - sampled reconstruction of size 256  256 as references .
a variable density cartesian sampling trajectory with 80% undersampling was employed as shown in figure 8(b ) .
the reconstructions from two data sets using dlmri , tbmdu and wtbmdu - l2 with p = 0.7 are displayed in figures 8(b ) , 8(c ) , and 8(d ) and figures 9(b ) , 9(c ) , and 9(d ) , respectively .
the corresponding psnr values of three methods are 31.35  db , 32.04  db , and 32.69  db for reconstructions in figure 8 and 30.58  db , 31.59  db , and 31.82  db for those in figure 9 .
some regions of error in these subplots have been illustrated with red arrows , indicating that the reconstruction using wtbmdu shows better fidelity to the reference image than that of dlmri and tbmdu .
the enlargements of the reconstruction results with these methods are shown in figures 8(f ) and 9(e ) .
it can be observed that visible aliasing artifacts along the phase encoding direction ( horizontal in the image plane ) in the dlmri reconstruction is more pronounced than those of tbmdu reconstruction , while there are almost no artifacts in that of wtbmdu with p = 0.7 .
close inspection of the zoomed - in images indicates that the weighted sparse regularization enhances the reconstruction quality . in summary
, the results from our work clearly demonstrate the advantages of introducing iterative reweighting scheme and alternating direction method in nonconvex optimization over conventional methods , for constrained image reconstruction from undersampled k - space data .
not only for the brain image containing large piecewise constant regions , as shown in figures 1 , 2 , and 3 , but also for the circle of willis and lumbar spine composed of many texture features as shown in figures 4 and 5 , the proposed method visually provides more pleasant results than existing methods . by means of the bregman iteration / al methodology ,
it is notable that the proposed method consistently supports superior results without any parameters tuned manually regardless of different sampling trajectories , varying sampling factors , and the existence of noise or not .
we note that , in some circumstances , the psnr value started to decrease a little bit after achieving the highest value .
the reason of this phenomenon is that the solution of problem ( 1 ) does not necessarily have higher psnr than the intermediate iterates . in real - time imaging applications where the acquisition and reconstruction speed is crucial , a numerical algorithm which converges fast at the beginning iterations
this work presents a new regularizer combining the nonconvex l 
 p pseudonorm and adaptive dictionary for mri reconstruction from undersampled k - space data .
based on the potential success of tbmdu in the al framework for cs - mri , the proposed wtbmdu method extends the plain tbmdu to handle nonconvex l 
 p regularization by weighting the sparse coding of the coefficients .
l 
 2 technique results in the wtbmdu - l1 and wtbmdu - l2 algorithms , respectively .
high - accuracy reconstructions are obtained from significantly undersampled k - space data by minimizing the adaptively sparsity - promoting regularization criterion subject to the data - consistency constraint . in particular , compared to the dlmri method that minimizes the l 
 0 regularization by the greedy algorithm of matching pursuit , the iteratively reweighted strategy combined with al framework yields much higher psnr values . on the other hand ,
compared to the plain tbmdu for l 
 1-norm minimization , the adaptively weighted method achieves 0.32.5  db psnr improvement and visually provides more pleasure results .
various experimental results demonstrate the superior performance of the method under a variety of sampling trajectories and k - space acceleration factors .
wtbmdu is very general since it represents an algorithmic framework that can be easily adapted to different reweighting strategies and nonconvex sparsity - inducing functions proposed in the previous literatures .
some extensions will be studied in future as follows : ( i ) for the reweighted l 
 1 strategy , using different nonconvex objective functions as conducted in [ 38 , 41 ] ; besides , the shrinkage can be generalized as in [ 47 , 48 ] ; ( ii ) for the reweighted l 
 2 strategy , trying different nonconvex objective functions as in .
additionally , the proposed scheme can be easily extended to incorporate other prior ( e.g. , spike and slab prior used in , where the regularization consists of the l 
 2-norm and l 
 0-norm terms ) .
it may be better for stabling the numerical process ; ( iii ) the current iteratively reweighted strategies can be combined with the homotopic technique described in [ 1416 ] , where the value of p slowly decreases from 1 to 0 .
another extension we will consider in future is to extend our proposed model to parallel mri reconstruction .
we note that a very recently published paper , under the work of ramani and fessler , addressed the regularized sense - reconstruction using al methods . in order to effectively solve the unconstrained sense - reconstruction consisting of tv / wavelets regularization term and data - fidelity term using the al formalism ,
they split not only the regularization term , but also the fourier encoding and spatial components in the data - fidelity term , by introducing auxiliary variables to decouple the data - domain components , and the regularization component , respectively .
we believe that the regularizer of the sparse representations of overlapping image patches in our work can be naturally incorporated into their work in a unified al framework .
 nonconvex optimization has shown that it needs substantially fewer measurements than l 
 1 minimization for exact recovery under fixed transform / overcomplete dictionary . in this work , 
 two efficient numerical algorithms which are unified by the method named weighted two - level bregman method with dictionary updating ( wtbmdu ) are proposed for solving lp optimization under the dictionary learning model and subjecting the fidelity to the partial measurements . by incorporating the iteratively reweighted norm into the two - level bregman iteration method with dictionary updating scheme ( tbmdu ) , the modified alternating direction method ( adm ) 
 solves the model of pursuing the approximated lp - norm penalty efficiently . 
 specifically , the algorithms converge after a relatively small number of iterations , under the formulation of iteratively reweighted l 
 1 and l 
 2 minimization . 
 experimental results on mr image simulations and real mr data , under a variety of sampling trajectories and acceleration factors , consistently demonstrate that the proposed method can efficiently reconstruct mr images from highly undersampled k - space data and presents advantages over the current state - of - the - art reconstruction approaches , in terms of higher psnr and lower hfen values . 
1. Introduction
2. Theory
3. Materials and Methods
4. Experiment Results
5. Conclusions and Discussions