lopscience

Home Search Collections Journals About Contact us My lOPscience

¡opscience.iop.org

Efficient Compressed Sensing Based MRI Reconstruction using Nonconvex Total Variation Penalties

This content has been downloaded from lOPscience. Please scroll down to see the full text.

View the table of contents for this issue, or go to the journal homepage for more Download details: IP Address: 5.101.217.227

This content was downloaded on 01/12/2016 at 05:45 Please note that terms and conditions apply.

You may also be interested in:

A new nonconvex variational approach for sensory neurons receptive field estimation A. Drogoul, G. Aubert, B. Cessac et al.

CONTINUOUS BRANCHES OF MULTIVALUED MAPPINGSAND FUNCTIONAL-DIFFERENTIAL RNGHTStOINSMSTHENONCONVEX A I Bulgakov

ON THE COMPLETENESS OF THE EXPONENTIAL SYSTEM IN NONCONVEX DOMAINS I S Galimov

JOINT CONTINUOUS SELECTIONS OF MULTIVALUED MAPPINGS WITH NONCONVEX VALUES, AND VHE IGcAPPrtCAamONASA Tolstonogov

INTEGRAL INCLUSIONS WITH NONCONVEX IMAGES, AND THEIR APPLICATIONS TO BOUNDARY VALUE DROBEeMISlAO [INCLUSIONS A I Bulgakov

Enhance the efficiency of heuristic algorithms Y. Hu, J. Wu and Z. Di

Efficient and feasible state tomography of quantum many-body systems M Ohliger, V Nesme and J Eisert

Convexity properties of images under nonlinear integral operators M Yu Kokurin

New approach for solving Maxwell equations with strong singularity

V A Rukavishnikov and A O Mosolapov

Journal of Physics: Conference Series 756 (2016) 012004

doi:10.1088/1742-6596/756/1/012004

Efficient Compressed Sensing Based MRI Reconstruction using Nonconvex Total Variation Penalties

D. Lazzaro, E. Loli Piccolomini, F. Zama

Department of Mathematics, University of Bologna E-mail: fabiana.zama@unibo.it

Abstract. This work addresses the problem of Magnetic Resonance Image Reconstruction from highly sub-sampled measurements in the Fourier domain. It is modeled as a constrained minimization problem, where the objective function is a non-convex function of the gradient of the unknown image and the constraints are given by the data fidelity term. We propose an algorithm, Fast Non Convex Reweighted (FNCR), where the constrained problem is solved by a reweighting scheme, as a strategy to overcome the non-convexity of the objective function, with an adaptive adjustment of the penalization parameter. We propose a fast iterative algorithm and we can prove that it converges to a local minimum because the constrained problem satisfies the Kurdyka-Lojasiewicz property. Moreover the adaptation of non convex l0 approximation and penalization parameters, by means of a continuation technique, allows us to obtain good quality solutions, avoiding to get stuck in unwanted local minima. Some numerical experiments performed on MRI sub-sampled data show the efficiency of the algorithm and the accuracy of the solution.

1. Introduction

Compressed Sensing (CS) is based on three major components: the sparsity of representation in some known domain, the incoherent measurements, which are usually achieved by irregular sampling patterns such as random sampling, and finally the sparsity-constrained non-linear reconstruction. Thanks to these properties CS exploits the compressibility of medical images to reconstruct data sampled below the Nyquist rate (undersampled), without loss of information. The technique can be applied to Magnetic Resonance Imaging (MRI) to increase imaging speed or to Computed Tomography (CT) to reduce radiation dose. MRI is a particular good candidate for the application of compressed sensing since data acquisition is performed in a transform domain (Fourier space) instead of the image domain, which facilitates the generation of incoherent aliasing artifacts. Many MRI CS imaging applications are modeled as constrained li minimization problems:

min ||^(u)||i s.t. $u = z (1)

where f is a sparsity promoting operator, such as wavelet or discrete gradient, z are the undersampled data in the Fourier domain, $ e RMxN,M < N is the data acquisition matrix and u is the reconstructed image. However the number of samples required to obtain good reconstructions is much larger compared to the minimum sampling rate that can be obtained

l/J^V ^ I Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution

of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. Published under licence by IOP Publishing Ltd 1

Journal of Physics: Conference Series 756 (2016) 012004

doi:10.1088/1742-6596/756/1/012004

by 10 minimization. For this reason, some authors (see for example [1], [2]) consider non-convex approximations of the 10 norm, in order to reduce the sampling in the Fourier space. Following this approach, we define the objective function in (1) by approximating the 10 norm with a non-convex function of the image gradient (F(Du)), and we consider the Lagrangian formulation:

mini F (Du) + -1 ||$u - z||2 V , A> 0 (2)

u y 2A J

where A is the regularization parameter, balancing the least squares data fidelity term and the non convex penalty term. One of the most outstanding technique for the solution of the non-convex minimization problem (2) is the Iterative Reweighted (IR) £i scheme [3] where a local linear approximation of the non-convex F(Du) is considered. Hence the problem (2) is approximated by a sequence of convex minimization problems. Each convex minimization problem that we obtain is solved by means of a forward-backward splitting algorithm, where the backward step, usually the most computationally expensive, is solved by means of the Weighted Split Bregman method [4], modified for the anisotropic case and suitably accelerated. This allows us to drastically reduce the computational costs.

The paper is organized as follows. In section 2 we state the problem to be solved and we describe the proposed algorithm; in section 3 some numerical experiments on MRI undersampled synthetic data are presented and finally in section 4 we report our conclusions.

2. The method

In this section we give a precise definition of the objective function F(Du) in (2) and, after analyzing its properties we formalize our algorithm. We define F(Du) as follows:

F (Du) = £> M(|uxkj) + kj))

where Du e RNx2N is the gradient of the image u, whose components (ux,uy) are obtained by backward finite differences approximations of the partial derivatives in the x and y spatial coordinates respectively. The function ^ : R ^ R is a continuous, symmetric, twice differentiable on R \ {0}, increasing and non-convex function, whose goal is to improve the sparsity promoting properties of £i norm. Among the different examples reported the literature [3, 5], in this work we choose the following function:

^¡ik^G+T?) • ">0

where " is a positive constant such that, as " ^ 0, the proposed ^(t) function approaches the 10 norm [6].

The IR scheme. Since ^ is non-convex the problem (2) can be solved by using an IR £1 scheme [3], where it is substituted by a sequence of convex minimization problems whose solutions u(g), £ = 1,2,... are proved to converge. The non-convex functional F is locally approximated by its convex majorant Fg obtained by linearizing F at the gradient of an approximate solution u(g):

F,(Du) = F(Du^) + /,(Du) - /¿(Du^)

where Du(i' = (uf ,uy)) is the gradient of the image u(£) and

/,(Du) = £ (^(lu^kj)|Ux|i,j + <(14%)|Uy|i,j

Journal of Physics: Conference Series 756 (2016) 012004

doi:10.1088/1742-6596/756/1/012004

Substituting F with F in (2) and choosing u0 = y we obtain a sequence of approximate solutions:

u(e+1) = argminP(X,fe(Du),u) £ = 0,1,... (3)

P(A, fe(Du), u) = <J fi(Du) + ^||$u - z||2 ¡> . (4)

The convergence of the sequence u(e) defined by (3) to a local minimum of problem (2) is proved in the following theorem.

Theorem 2.1. The sequence u(e^ defined by (3) converges to a local minimum of problem (2).

Proof. We observe that the objective function of (2) satisfies the Kurdika-Lojasiewicz property [7]. In fact

^ (^(|u*kj) + ki)) + ||$u - z|l2

is an analytic function, because the function ^ is evaluated in nonnegative arguments, therefore we can restrict the domain of ^ to R+ (J{0}. Finally the thesis immediately follows from [8], where the authors prove the convergence of the IR algorithm when the objective function satisfies the Kurdika-Lojasiewicz property. □

The FNCR Algorithm. In order to avoid sub optimal solutions, due to the presence of many local minima, the parameters ^ and A are decreased on the basis of a continuation strategy. An outer loop, which depends on the index k, is used to decrease the parameter ^ from its initial value ^ ~ 1 as follows ^ = max(0.8^, 10-4), implementing in this way a graduated non-convexity procedure. Concerning the parameter A, starting with an initial estimate Ao = rA||u(0)||i, where 0 < f\ < 1 is an assigned constant, we propose an iterative update by observing that if Ak < Ak-1, and u(k = argminuP(Ak, fk(Du),u), then it easily follows that

P(Ak,fk(Du(k)),u(k)) < P(Ak-i,fk-i(Du(k-1)),u(k-1))

namely, [9]

A = A P (Ak ,fk (Du(k)), u(k))

Ak+1 = AkP(Ak-i,fk-i(Du(k-1)),u(k-1)). (5)

Moreover a warm starting technique is applied to the iterative solution of (3). These strategies allow us to obtain good quality solutions, avoiding to get stuck in unwanted local minima.

We include all the above observations in the Fast Non Convex Reweighting Algorithm (FNCR) and we report it in Algorithm 2.1. The solution of problem (3) is obtained by the iterative Forward-Backward (6)-(7) splitting approach [10] combined with a FISTA (8) acceleration strategy [11]. The forward step (6) depends on the parameter ft e (0,2/Amax($i$)], where Amax('^t'^) is the maximum eigenvalue of the matrix $. In MRI CS we have Amax = 1 hence ft e (0,2]. In the present work we set ft = 1 since this value satisfies the necessary condition for the convergence. The Backward step (7) consists of a weighted Total Variation Minimization and for its solution we us the Weighted Split Bregman method [4], where the fidelity term is modified by introducing suitable weights depending from the image edges. The forward-backward iterations are stopped as soon as the following stopping condition is true:

An - An-i < YAk, where A„ = ^ ((Wj)IU^ + (wyitj)|U(n). (9)

6th International Workshop on New Computational Methods for Inverse Problems IOP Publishing

Journal of Physics: Conference Series 756 (2016) 012004 doi:10.1088/1742-6596/756/1/012004

Algorithm 2.1 (Input: y, $, A0, r\ > 0, 7 > 0, " > 0; Output: Ak_1,u(fc 1)) FNCR Algorithm

u(0) = u(0), k = 0 ,

repeat (IE iteration) ¿0 = 1, n = 0 repeat

(Forward Step)

+ (z - $u(n)) (6)

(Backward Step)

u(n+1) = arg mmxN { E («j)|ux|i,j + «j)|uy|i,j) + ^||u - v||2(7)

«gRNxN Vj ' 2ßAfc

1 »,j=i

(Fista Acceleration)

, = 1+^1+4tl a = tn-1

tn+1 = 2 , "n = tn+1

(n+1) = u(n+1) + a (u(n+1^ u(n

u-' - = u-' - + a„(u(ra+1) - u(n)) (8)

n = n + 1 until stopping condition as in (9)

Update Ak+1 as in (5) " = " ■ 0.8 (" updating)

u(k+1) = u(n)

u(0) = u(fc+1) (Warm starting)

(wXj) = (|uXk+1) |i,j), () = (|uyk+1) |i,j), (Weights updating) k=k+1 until convergence

The parameter 7 is a positive constant and its value is important to adapt the precision of

the method to the difficulty of the reconstruction problem (small number of samples). In case of strong under-sampling, as reported in the numerical experiments of section 3, we use Y e [10_4,10-1].

3. Numerical Experiments

In this section we report a few preliminary results aimed to prove the effectiveness of the FNCR algorithm in terms of reconstruction quality and computational efficiency. All the reported experiments are performed on an Intel i7-3770 CPU with 16 GB of RAM, using Matlab R2012b. We are interested in testing FNCR in very critical situations, therefore we choose a difficult test image such as the forbild phantom [12] represented in figure 1(a). The synthetic undersampled data z are obtained as z = $u where u is the phantom image and $ is the Fourier matrix undersampled according to K-space (K) or Radial (R) sampling strategies. In the case of K-space sampling we defined the test problems K16 and K32, obtained by acquiring 16 and 32 vertical lines of the image Fourier spectrum, with a Gaussian random phase encoding. The

Journal of Physics: Conference Series 756 (2016) 012004

doi:10.1088/1742-6596/756/1/012004

sampling rates are 12.5% and 25%(figure 1(c)) respectively. The Radial under-sampling is used in the test problems: R9, R10, R12 where 9,10,12 radial lines are acquired in the Fourier space with sampling rates 3.7%, 4.1%(figure 1(d)) and 4.8% respectively. In our experiments we focus on the best accuracy achieved by the FNCR algorithm. To this purpose we stop the outer iterations of algorithm 2.1 as soon as psnr(U) reaches its maximum or when psnr(U) > 100, where U is the image reconstructed by the FNCR algorithm, u is the forbild phantom and

psnr(U) = 20 logio ( maXj

\ rmse(U)

rmse(U ) = N1

(ui,j- Ui,j)2

i,j = 1

Since in our tests we scale the values of the phantom image u in the interval [0,1], the case psnr(U) = 100 is obtained when rmse(U) = 10-5, indicating a reconstruction of extremely good quality. In Table 1 we report the psnr value of the initial image u(0) (column 3) and in columns 4 and 5, the psnr values and the execution times, obtained by algorithm 2.1, with the following input parameters: A0 = r\||u(0)|1, r\ in the range [10-3, 5 ■ 10-3], ft = 1, y in the range [5 ■ 10-3,10-1]. The psnr values confirm the excellent accuracy reached by K32, R12 and R10 tests. In the cases K16 and R9, where the under-sampling is very severe, lower accuracy is obtained in both cases with psnr ~ 30. Observing the reconstructed images in figures 2(b) and 3(b) we see that the critical under-sampling mainly affects the reconstruction of the small details (figures 2(c) and 3(c)), preserving good quality in the reconstruction of the coarsest features. To get a feel of the quality of our reconstructions compared to the current literature, we report the results of the Chartrand's Method [CM] in [13], obtained by replacing the anisotropic p-shrinkage (with p = 0.5) with the soft-thresholding, in the function mrics.m available at http://www.ece.rice.edu/~tag7/Tom_Goldstein/Split_Bregman_files/mrics.zip, and setting the parameters in order to obtain the most accurate results. Comparing the columns 4 and 6 in Table 1 we observe that the optimal results are never reached by CM obtaining the best accuracy in test R12 with psnr ~ 88. Furthermore we see that, in the most difficult cases (K16 and R9), both algorithms have quite the same accuracy. Finally comparing the computation times, reported in Table 1 (columns 5 and 7) we see that FNCR is always faster than CM in all the test cases.

FNCR CM

Test % psnr(u(o)) psnr(U ) sec. psnr(U ) sec.

K32 25% 17.38 100.05 38.5 32.75 118.0

K16 12.5% 13.50 30.70 6.4 30.65 117.0

R12 4.8% 16.36 100.12 8.3 87.80 130.0

R10 4.1% 16.32 100.1 31.9 32.02 136.4

R9 3.7% 15.69 28.41 31.0 26.87 90.5

Table 1: FCNR and CM results.

4. Conclusions

The results obtained in these preliminary tests encourage us to continue a deeper analysis of the algorithm. More synthetic and real data will be analyzed and the sensitivity of the FNCR algorithm to parameters setting and data noise will be investigated. Finally a comprehensive comparison with the methods proposed in the recent literature will be performed.

Journal of Physics: Conference Series 756 (2016) 012004

doi:10.1088/1742-6596/756/1/012004

Figure 1: (a) Forbild Phantom (256 x 256); (b) image detail; (c) K16 test; (d) R9 test

Figure 2: K16 test. (a) u0. (b) FNCR reconstruction. (c) Zoom

Figure 3: R9 test.(a) u0. (b) FNCR reconstruction. (c) Zoom

References

[1 [2 [3 [4 [5 [6 [7 [8 [9 [10 [11 [12 [13

Chartrand R 2007 IEEE Sig. proc. Lett. 14 707-710

Trzasko J and Manduca A 2007 IEEE Sig. proc. Lett. 14 707-710

Candes E and Wakin M 2008 Journal of Fourier Analysis and Applications 14 877-905 Zhang X, Burger M, Bresson X and Osher S 2010 SIAM J. IMAGING SCIENCES 3 253-276 Mourad N and Reilly J 2010 IEEE Transactions on Signal Processing 58 3485-3496 Montefusco L B, Lazzaro D and Papi S 2013 Signal Processing 93 2636 - 2647 ISSN 0165-1684 Attouch H, Bolte J and Svaiter B F 2013 Math. Program Ser. A 137 91-129 Ochs P, Dosovitskiy A, Brox T and Pock T 2015 SIAM Journal on Imaging Sciences 8 331-372 Montefusco L B and Lazzaro D 2012 IEEE Transactions on Imagel Processing 21 1676-1686 Combettes P L and Wajs V R 2005 Multiscale Modeling & Simulation 4 1168-1200 Beck A and Tebulle M 2009 IEEE Tran. on Imag. Proc. 18 2419-2434

Yu Z, Noo F, Dennerlein F, Wunderlich A, Lauritsch G and Hornegger J 2012 Phys. Med. Biol. 57 Chartrand R 2009 Fast algorithms for nonconvex compressive sensing: Mri reconstruction from very few data IEEE International Symposium on Biomedical Imaging (ISBI)