magnetic resonance imaging ( mri ) is one of the major imaging modalities in use today . compared to computed tomography ( ct ) , mri has advantages in imaging soft tissues .
however , its relatively long imaging time remains a great challenge for clinical application , often limiting its application .
significant efforts have focused on faster data collection as well as reducing the amount of data required without degrading image quality .
for example , parallel imaging [ 13 ] exploits redundancy in k - space by introducing multiple receiver channels , mitigating the aliasing artifacts caused by a reduced sampling rate .
recently , compressed sensing based mri ( cs - mri ) allows high quality reconstruction from undersampled data by enforcing the pseudo - sparsity of images in a predefined basis or dictionary , such as the traditional two - dimensional ( 2d ) separable wavelet transform or total variation .
the discrete wavelet transform ( dwt ) , for example , has three major disadvantages : shift sensitivity , poor directionality , and lack of phase information [ 7 , 8 ] . for these reasons , traditional dwts fail to capture regularities of contours , since they are not able to sparsely represent one - dimensional singularities of 2d signals .
therefore , improvements can be obtained by mitigating some of these disadvantages simultaneously . in this paper
, we propose an improved compressed sensing method for mr imaging by utilizing the double - density dual - tree dwt .
the authors in   used dual - tree complex wavelet transform ( dt - cwt ) as a sparsifying transform , which only has wavelets oriented in six directions .
but as natural images exhibit smooth regions that are punctuated with edges at several orientations , dual - tree complex wavelet transform may fail to sparsely represent the geometric regularity along the singularities , which require higher directional selectivity .
but they can only sparsely represent the smooth contours but not the points in images [ 13 , 14 ] . in this paper
, we propose one possible solution by using a newly developed multiresolution tool , double - density dual - tree transform , which may provide sufficient sparse representation for mr images with different features .
total variation is also exploited as a penalty in the reconstruction formulation to suppress noise .
note that in , the authors applied their method to radial trajectories . in this paper ,
its variant on cartesian sampling is used for comparison , namely , cs dt - cwt . to differentiate between the original compressed sensing based mri algorithm   and the improved version
, we will denote the original algorithm in   as cs and our proposed method as ics ( improved compressed sensing mri ) .
, we briefly review principles of mri and then discuss the design of our proposed algorithm . in section 3 , we will present the experimental results of our algorithm in comparison with some other algorithms . finally , a brief conclusion will be drawn .
  mri signal is generated by the proton in hydrogen atoms , the main component of the human body .
since protons are charged particles , when a human body is placed in an strong static magnetic field b0 , protons will align themselves with the magnetic field , yielding a net magnetic moment precessing around b0 .
the frequency of larmor precession is proportional to the applied magnetic field strength as defined by 
 ( 1)f=b0 , 
 where  is a constant ( 42.57  mhz / t ) .
if the frequency of the applied pulse is equal to the larmor frequency , the net magnetic moment will tilt away .
once the rf signal is removed , the protons realign themselves such that the net magnetic moment is again aligned around b0 .
the protons return to equilibrium by emitting rf signal , which is then captured by a conductive field .
this measurement is reconstructed to obtain gray - scale mr images . to produce a 3d image , a gradient magnetic field , gz
, is added to b0 so that the larmor frequency changes linearly in the axial direction , z. hence , an axial slice can be selected by choosing a specific larmor frequency of that slice .
additionally , two gradient magnetic fields , gx and gy , are applied causing the resonant frequencies of the protons to vary according to their positions in the x - y plane . as a result
if the signal is fully sampled at the nyquist rate , a 2d inverse fourier transform is then used to transform the encoded image to the spatial domain .
consider the following : 
 ( 2)x = f1y , 
 where y is the measurements from scanner , which is also called k - space data , f is the fourier transform matrix , and x is the desired mr image .
but in the real world , downsampling may be needed for some applications , such as to fit the scans into one - breath hold or to enable real time - imaging .
the 1d double - density discrete wavelet transform ( dd - dwt ) is based on one scaling function ( i.e. , low pass ) and two different wavelets ( i.e. , high pass ) where one wavelet is a half - sample shift of the other .
the 2d transform applies the 1d transform alternately to the rows and column , giving nine subbands , where one is a low pass subband and the remaining eight subbands become eight wavelet filters .
thus , the 2d complex double - density dual - tree discrete wavelet transform ( cdddt - dwt ) is implemented by using four 2d dd - dwt in parallel with different filter banks for rows and columns separately .  
lf and hf are the filter banks of the first level decomposition , which represent one scaling filter and eight wavelet filters .
l and h are the filter banks for the second and remaining levels of decomposition .
therefore , four low pass subbands ( l11l14 ) and 32 high pass subbands ( h11h14 ) are produced after one level of the transform . as a result ,
32 wavelets are created by the sum and difference on each pair of subbands . 

the figures are obtained by setting all the wavelet coefficients to zero , for the exception of one wavelet coefficient in each of the high pass subbands of one level .
therefore , if the transform has more directional wavelets , then fewer coefficients are needed to represent a given geometric object .  
obviously , 2d wavelet transform can resolve only three spatial - domain feature orientations : vertical , horizontal , and diagonal .
in addition , the third wavelet does not have a dominant orientation , which is the main cause of artifact ( checker board pattern ) . therefore
, traditional 2d separable wavelet fails to sparsely represent geometric structures , such as edges . in figure 2(b ) , six orientations ( 15 , 45 , 75 ) are obtained by 2d dt - cwt . by contrast
, complex double - density dual - tree wavelet transform has 32 wavelets oriented in 16 different angles ( see figure 2(c ) ) .
the 16 wavelets shown in the first ( second ) row can be interpreted as the real ( imaginary ) parts of the set of 16 complex wavelets and the third row are the magnitudes of the 16 complex wavelets
additionally , all of the wavelets are free of checker board pattern . as a result
, we can get more accurate representation for the geometric regions .   in figure 3
the first image group only contains a curve purposely designed to demonstrate the improved directionality property of cdddt - dwt .
we reconstruct the test images using the first stage of transform coefficients . as we can see from figures 3(b ) and 3(f ) , there are blocking artifacts in dwt reconstructions since it can only accurately represent vertical and horizontal lines . because dt - cwt has more directional wavelets , the reconstructed curve looks smoother and artifacts are reduced
. but six orientations ( figure 2(b ) ) are not sufficient to accurately represent this curve as it contains all directions .
this is because it has more wavelets which are strongly oriented at 16 different angles .
one may also conclude that 16 orientations are sufficient to represent any geometric object with high precision .
compared to the curve , the points are better reconstructed by all the three transforms . as stated in the introduction section , wavelet transforms are good at capturing point singularities . from this example , we can see that cdddt - dwt can sparsely represent both contours and points , indicating its superiority for these tasks .
the problem of undersampling k - space data actually leads to an underdetermined system of linear equations .
one way to improve performance is to incorporate a prior knowledge into the reconstruction process which is based on the idea of sparsity in compressed sensing ( cs ) [ 2123 ] .
the essence of cs is that a signal , which in our case is the image x , can be completely reconstructed with a high probability with far less samples than required by conventional nyquist - shannon sampling theorem , if the image has a sparse / compressible representation in a transform domain  , such that most entries of the vector x are zero or close to zero .
the entire process consists of three steps : encoding , sensing , and decoding . in the first step , the object image x of size
n is encoded into a smaller vector y = fux of a size m ( m < n ) by the system matrix .
fu denotes the fourier matrix associated with some undersampled trajectory while y is the corresponding undersample k - space data . directly solving
the image is approximately sparse in cdddt - dwt , that is ,  = x , has few elements with relatively large magnitudes , a solution is possible .
then the second step is obtaining the undersampled k - space data y from the imaging system . incorporating the sparsity prior knowledge into the process of image reconstruction
, the third step is to recover  ( and thus x ) by solving the following constrained optimization problem : 
 ( 3)argminx||x||1  subject  to  ||fuxy||2< , 
 where  is a parameter controlling the data consistency .
it has been mathematically proven that if the image has only k entries with relatively large magnitudes , the order of klog(n ) measurements is sufficient to accurately reconstruct x via the 1   norm minimization procedure with high probability . as we have noticed from figure 3 ,
therefore , in our proposed algorithm , we include the total variation ( tv ) as a penalty because it was shown that it is efficient in suppressing the noise in the reconstructed image .
the constrained problem in ( 3 ) can also be converted into an unconstrained problem , giving rise to our proposed ics model : 
 ( 4)arg  minx1||x||1+2||x||tv+||fuxy||22 , 
 where two regularization factors 1 and 2 are introduced to leverage the cost function 's emphasis on the transform 1   penalty , tv penalty , and the data fidelity term .
the selection of regularization factor has been an interesting area of research in the field of regularized iterative methods [ 25 , 26 ] .
a large 2 tends to suppress image gradient and make the reconstructed image look smooth , losing point - like features . in our study
, we chose the optimized regularization parameters 1 and 2 for all methods for a fair comparison .
the tv term of an image in this work is defined as follows : 
 ( 5)||x||tv=|x|dx . in a discrete version , ( 5 ) becomes 
 ( 6)||x||tv=i , j(xi+1,jxi , j)2+(xi , j+1xi , j)2 .  
to speed up the implementation , we exploit a fast implementation of cdddt - dwt .
since ( 4 ) poses an unconstrained convex optimization problem , we propose solving it using a nonlinear conjugate gradient descent algorithm that is similar to .
it has been shown in   that the iterative algorithm in our study has better performance than the algorithm in .
the conjugate gradient requires the computation of j(x ) which is 
 ( 7)j(x)=1||x||1+2||x||tv+2fu(fuxy ) .
as the 1 norm and total variation term ( 5 ) is the sum of absolute values , the absolute values , however , are not smooth functions , and as a result ( 7 ) is not well defined . in , lustig et al .
approximated the absolute value with a smooth function , |x|xx+ , where  is a positive smoothing parameter .
in this section , we report our experiments to evaluate and validate the proposed algorithm .
there are five sets of experiments . in the first two experiments , numerical phantoms and simulated k - space data
the third and fourth experiments used real data collected from real scanners . in the fifth experiment
, we manually add noise to the k - space data of the fourth data set .
the first phantom that we consider is the discrete shepp - logan , which contains geometrical structure and directional - oriented curves .
the second phantom is purposely designed to be nonsparse under total variation domain and contains features difficult to reproduce with partial fourier sampling . in such way
the final dataset was a coronal section of a brain obtained from a t1-weighted brain scan . to have a clear differentiation between these datasets , we name these four datasets as shepp - logan , phantom , axial brain , and coronal brain , respectively , as shown in figure 4 .
it should be pointed out that the proposed method also works with non - cartesian sampling pattern , such as spiral and radial trajectories , although these are not shown .
reconstructions were performed at different sampling rates : 0.1 , 0.15 , 0.2 , 0.25 , 0.3 , 0.35 , and 0.4 .
a sampling pattern with sampling rate 0.2 is shown in figure 5 where white bar means that the data is sampled and black means otherwise ( i.e. , not sampled ) .
the reconstructed data was quantitatively evaluated in terms of peak signal - to - noise rate ( psnr ) and structural similarity ( ssim ) index .
psnr measures the difference between the reconstructed image x^ and original image x , which is defined by 
 ( 8)psnr=20  log10(maxmse ) , 
 where mse=(1/n)i=0n(xi - x^i)2 and max is the maximum possible pixel value of the image .
the structural similarity ( ssim ) index is highly effective for measuring the structural similarity between two images .
suppose  and t are local image patches taken from the same location of two images that are being compared .
the local ssim index measures three similarities of the image patches : the similarity of luminances l( , t ) , the similarity of contrasts c( , t ) , and similarity of structures s( , t ) .
local ssim is defined as 
 ( 9)s(,t)=l(,t)c(,t)s(,t )  = ( 2t+c1p2+t2+c1)(2t+c2p2+t2+c2)(2t+c3t+c3 ) , 
 where  and t are local means ,  and t are local standard deviations , and t is cross - correlation after removing their means .
the higher the value of ssim , the higher image quality is delivered . in this paper , three methods will be compared under identical conditions .
the first method uses the discrete wavelet transform ( shown as cs ) ; the second method uses dual - tree complex wavelet transform ( shown as cs dt - cwt ) ; the proposed method uses complex double - density dual - tree wavelet transform ( cdddt - dwt , denoted by ics for the rest of the paper ) .
reconstruction was done under the same conditions such as iterative algorithm and sampling pattern as the accuracy depends on the selection of optimum regularization parameters .
we have used the last dataset , coronal brain , as an example to show the methodology of determining the optimal parameters . for each algorithm
then we set 2 = 0.001 and searched the optimal value for 1 that gives the highest psnr , as shown in figure 6(b ) .
thus , we used this recurring process to determine the optimum values of 1 and 2 .
figures 7(b ) , 7(c ) , and 7(d ) are reconstructions by cs method , cs dt - cwt , and our proposed method ics , respectively .
the weakness of the traditional separable wavelet transform in representing two - dimensional singularities , for example , curves and edges , is visualized with the reconstructed shepp - logan images . as we have seen that separable wavelet transform can resolve only three spatial - domain feature orientations .
therefore , the ellipses are not sparsely represented , causing substantial artifacts in the reconstructed images ( figure 7(b ) ) .
for example , two small ellipses around center ( as indicated by white arrows ) are faded away from the cs reconstruction . to see it clearly , detailed zoom - ins
by contrast , the proposed method reconstructs the image with higher visual image quality . from figures 7(e)to 7(h )
, we can see that , with more directional wavelets , the ellipses are better reconstructed .
the reconstruction by complex double - density dual - tree wavelet transform has the best visual image quality .
additionally , the two small ellipses could be clearly recognized ; see figure 7(h ) .  
similar results can be obtained from the second phantom experiment . from figure 8(b ) , we can see significant aliasing artifacts in the cs reconstruction .
these artifacts are caused by the use of nonideal low pass and high pass filters in the traditional separable wavelet transform .
additionally , from the detailed zoom - in ( figure 8(f ) ) , we can see that cs can not reconstruct those closely placed lines separately .
one possible reason for this phenomenon may be that the wavelets of dw - cwt and cdddt - dwt are finer than those of dwt see figure 1 .
this is evident that our method leads to a better reconstruction with higher spatial resolution . in order to see how the results vary with the sampling rates ,
experiments were also performed at sampling rates of 0.1 , 0.15 , 0.25 , 0.3 , 0.35 , and 0.4 .
the plots of psnr and ssim versus different sampling rates for the phantom simulations are shown in figure 9 . at sampling rate 0.25 ,
the psnr of the shepp - logan image by our proposed method is 41.5  db , and structural similarity is well above 0.93 .
there is no visual difference between the reference and reconstructed images using the proposed method . on the other hand ,
reducing the sampling rate below 0.2 , reconstructed images are obviously blurred and their quality is not acceptable .
one possible reason for the failure of compressed sensing methods is that the initial image is too bad at such low sampling rates .
figures 10 and 11 show the reconstructed images where ( a ) , ( b ) , ( c ) , and ( d ) are original images , cs reconstructions , cs dt - cwt reconstructions , and our ics reconstructions , respectively .
note that although both cs reconstructions are able to reduce the undersampling artifacts significantly , the images obtained by cs contain significant blocking artifacts along directional edges . to emphasize this point , detail views of the reconstructed images
this is because the diagonal wavelet of traditional wavelet transform does not have a dominant orientation see figure 1(a ) . as a result
, it can not represent the curves efficiently , leading to considerable artifacts along edges .
dual - tree complex wavelet transform can resolve six orientations , both of which have a dominant orientation , and thus the image quality is greatly improved by cs dt - cwt . in contrast , since complex double - density dual - tree wavelet has improved directional selectivity and finer wavelets than dwt and dt - cwt , the proposed method preserves more edges and small structures . therefore , the images obtained using our proposed method are much closer to the original images .   figure 10(h ) and
figure 11(h ) indicate that edges are sharper in the reconstructed image using our proposed method .
the evaluation matrices ( psnrs and ssims ) of the reconstructed images are also plotted against the sampling rates in figure 12 .
figure 12 also indicates that psnrs and ssims of reconstructions by our proposed method in all cases are higher than those of other methods .
one should also note that the original cs method can produce an accurate reconstruction given a sufficient sampling rate .
for instance , the results in   showed that cs can produce high quality brain image without artifacts at sampling rate 0.42 ( 2.4 acceleration ) .
note that the psnr of cs method is well above 35  db and ssim is as high as 0.9 at sampling rate 0.4 . but at low sampling rates , reconstructed image is not satisfactory . as a result
our proposed complex double - density dual - tree wavelet is able to capture those geometric features that are not captured by other transforms , producing a higher quality of image .
the results of this test confirm that our proposed method outperforms the original cs method in preserving edges and suppressing undersampled artifacts .  
finally , we reconstructed the coronal brain image using imperfect k - space data to test the robustness of the proposed method .
additive gaussian white noise e of relative magnitude ||e||2/||fxture||2 = 0.05 was purposely added to the k - space data .
, we can see that the parameter 2 is bigger when noisy data is used for reconstruction for all the methods .
but if 2 is too large , the reconstructed image will be oversmooth , losing low contrast or detailed information .
although the image quality is reduced , the proposed method can still maintain the best results . from the various results presented in this section , proposed method delivers higher psnr and ssim than other methods .
however , the computation time of cs - based methods including the proposed algorithm is about 70 seconds .
development of faster algorithms for solving problem ( 4 ) will be pursued in the future work .
in this paper , an improved compressed sensing mri method is proposed . the directionality property of complex double - density dual - tree wavelet transform is investigated . considering the edge features and quality , we employ objective measurements to evaluate the performance of our approach .
simulation results on phantom and in vivo data demonstrate that the proposed method can better reconstruct the edges and reduce undersampled artifacts than traditional cs - mri method does . in our implementation
, we use the nonlinear conjugate gradient iterative method to solve the unconstrained optimization problem .
further effort is needed to utilize a more advanced iterative method to improve reconstructed precision and reduce computational time for real time - imaging .
extension of this work with the nonconvex optimization is being considered to further improve the reconstruction .
 undersampling k - space data is an efficient way to speed up the magnetic resonance imaging ( mri ) process . 
 as a newly developed mathematical framework of signal sampling and recovery , compressed sensing ( cs ) allows signal acquisition using fewer samples than what is specified by nyquist - shannon sampling theorem whenever the signal is sparse . as a result , cs has great potential in reducing data acquisition time in mri . in traditional compressed sensing mri methods , 
 an image is reconstructed by enforcing its sparse representation with respect to a basis , usually wavelet transform or total variation . in this paper 
 , we propose an improved compressed sensing - based reconstruction method using the complex double - density dual - tree discrete wavelet transform . 
 our experiments demonstrate that this method can reduce aliasing artifacts and achieve higher peak signal - to - noise ratio ( psnr ) and structural similarity ( ssim ) index . 
1. Introduction
2. Theory and Method
3. Numerical and Experimental Results
4. Conclusion