We consider the reconstruction of MRI images by minimizing regularized cost-functions. To accelerate the computation of the estimate, two forms of half-quadratic regularization, multiplicative and additive, are often used. In Nikolova and Ng (2002), we have compared both theoretically and experimentally the efficiency of these two forms using one-dimensional signals. The goal of this paper is to compare experimentally the efficiency of these two forms using MRI image reconstruction. We find that using the additive form is more computationally effective than using the multiplicative form.published_or_final_versio

Deriche, R

Kornprobst, P

Ng, KP

Nikolova, M

HKU Scholars Hub

Title Half-quadratic regularization for MRI image restorationAuthor(s) Deriche, R; Kornprobst, P; Nikolova, M; Ng, KPCitationIEEE International Conference on Acoustics, Speech and SignalProcessing Proceedings, Hong Kong, China, 6-10 April 2003, v.6, p. 585-588Issued Date 2003URL http://hdl.handle.net/10722/46609Rights©2003 IEEE. Personal use of this material is permitted. However,permission to reprint/republish this material for advertising orpromotional purposes or for creating new collective works forresale or redistribution to servers or lists, or to reuse anycopyrighted component of this work in other works must beobtained from the IEEE.HALF-QUADRATIC REGULARIZATION FOR MRI IMAGE RESTORATION Rachid Deriche, Pierre Kornprobst INRIA 2004 Rte des Lucioles BP 93, 06560 Sophia-Antipolis, France ABSTRACT We consider the reconstruction of MRI images by minimiz- ing regularized cost-functions. To accelerate the computa- tion of the estimate, two forms of half-quadratic regulariza- tion, multiplicative and additive, are often used. In [13], we have compared both theoretically and experimentally the ef- ficiency of these two forms using one-dimensional signals. The goal of this paper is to compare experimentally the effi- ciency of these two forms using MRI image reconstruction. We find that using the additive form is more computation- ally effective than using the multiplicative form. 1. INTRODUCTION We address image reconstruction where a sought image ? E RP is estimated from degraded data y E Rq by minimizing a cost function J : Rp --f R which combines a quadratic data- fidelity term and a regularization term Q via a parameter i = min J ( z ) ;  where J ( z )  = 11Az-g112+BQ(z). (1) We shall assume that the observation operator A E RqXp is known. We focus on regularization term p > 0: ZERO of the form ,= 1 where $ : R --f R is a potential function and yT. for i = 1,. . . , r ,  are linear operators. Typically, {gTz} are first or second-order differences between neighboring pixels. If G THE WORK HAS BEEN SUPPORTED BY A GRANT FROM THE PROCORE - FRANCEIHONG KONG JOINT RESEARCH SCHEME SPONSORED BY THE RESEARCH GRANTS COUNCIL OF HONG KONG AND THE CONSULATE GENERAL OF FRANCE IN HONG KONG (REFERENCE NO. F ~ H K  30/01C. Mila Nikolova DCpartment Traitement du Signal et des Images, ENST 46, rue Barrault, 75634 Paris CCdex 13, France Michael N g  Department of Mathematics The University of Hang Kong Pokfulam Road, Hang Kong is the r x p matrix whose ith row is gT. for i = 1 , .  . . , r ,  a basic requirement is ker(ATA) fl ker(GTG) = {O}. We suppose that $ is smooth and convex, and edge-preserving, i.e. $( t )  < tZ  for ltl + 00. Such functions can be found in [2 ,  3, 1, 51, e.g., Huber potential function: (3) if It1 5 U ,  Cost-functions of this form are popular in various inverse problems such as denoising, deblurring, seismic imaging, tomography. However, the resultant minimizers ? are non-linear with respect to data y their computation is costly, especially when A has many non-zero entries. In order to cope with nu- merical slowness, halfquadratic (HQj reformulation of J has been pioneered, using two different ways, in [7] and [SI. The idea is to construct an augmented cost function J' : Rp x R" + R which involves an auxiliary variable s E R', and two new functions, Q : R x R --t R ,  where & ( . , s i )  isquadraticVsiER,Bnd1C. : R + R ,  J ' ( z , s )  = I I .hz-~l l~ +PcQ(~r~,si)+Bc~i)(~i), (4) so that $(t) = min {Q(t,  s) + $J(s)} , Vt  E R. ( 5 )  By (5 ) .  the global minimizer (?,2) of J' yields the solu- tion initially defined in (1 j, since J ( z )  = minSER. J ( z ,  s), Vz E Rp. In [7], Geman & Reynolds first considered a quadratic term Q of the multiplicative form, Later, Geman & Yang [XI proposed an additiveform for Q: (7) In both cases (6) and (7j, the dual function $, which ensures (5j, is obtained using the theory of convex conjugacy [9]. i=l i=l SER Q(t , s )  = t 2 s ,  for t E R, s E R+. (6) Q(t , s )  = ( t  - s)~, for t E R, s E R. 0-7803-7663-3/03/$17.00 02003 IEEE VI - 5x5 ICASSP 2003 The augmented cost-function 9 is minimized using an alternating minimization scheme.  Let the solution obtained at iteration ( k  - 1) read ( d - ' ) ,  dS-l)). At the next itera- tion k we calculate s ( k )  such that J ( Z ( ~ - ' ) , S ( ~ ) )  < - 9 ( d k - ' ) , s ) ,  Vs E R", dk)  such that Vx E Rp. These minimizations give rise to tn'o minimizer mappings,  x 4 [u(gTx),  . . . ,u(gTz)IT with U : R + R, and s + x ( s )  with x : R' 3 RP. The alternate minimization thus .7(dk), dk)) 5 J ( z ,  .dk)) ,  #) = x(s (" ) .  (9) These ideas has been pursued and deepened by many authors [3,6, 1, 5,  11, 101. Although the intuition that HQ regularization does indeed increase the speed of the mini- mization of regularized cost-functions of the form ( l ) ,  this critical question has never been considered in a theoreti- cal way. In [13], the performance of both formulations (6) and (7) has been compared using one-dimensional signals. The goal of this paper is to compare experimentally the effi- ciency of these two forms using MRI image reconstruction. We find that using the additive form is more computation- ally effective than using the multiplicative form. 2. S O M E  FACTS ABOUT HQ REGULARIZATION 2.1. Multiplicative form We consider potential functions $ such that t + $(&) is concave on R+, limtbo $'(t)/ t  <CO, t + $( t )  is convex on R, limt,, $(t)/ t2 = 0, 4 is twice differentiableon R, $(t) =$(-t), VtER. Then the expressions below are equivalent [7,S, 101: (lo) (1 1) $(t)  = inf,€R {st' + $(s)} , Notice that 21, is convex and $(s) = +CO f o r s  < 0; hence the infimum in ( I  1) can be considered only for s > 0. The resultant augmented cost-function 3 is defined on RY x R; and reads J ( s , s )  = I IAx-~11 '+P(Gz)~diag(s )  G z + p X $ ( s i ) ,  (12) where diag(s) is a diagonal matrix whose diagonal elements are si, for i = 1;. . . , r .  The function U ,  as given in (81, reads 15. 101 = suPt€R { $ ( t )  - s t ' }  ' i= 1 . .  ~ 1 3 )  where clearly u(t)  2 0, Vt E R .  The minimizer mappin 2, introduced in (9), satisfies D l J ( x ( s ) ,  s) = 0, Vs ER', and reads: (14) X ( S )  = (H(s) ) - '  ATy, where H ( s )  = A T A  + BGTdiag(s)G. 2.2. Additive form This form is considered under the condition that the func- is convex, continuous and finite for every t E R. Then the following expressions are equivalent: $( t )  = i&x {$(SI + (t - s)'/2}, (16) $(s) = SUPtCR {$ ( t )  - ( t  - s ) 2 / 2 } .  The condition ( I S )  implies that $'( t - )  2 @(t+), for any t E R. Whenever 4 is convex, it implies that $ is differen- tiable. The augmented cost-function now reads P J ( z ,  S) = IIAz-?/I12 + -IIGz -SI/' + B  $(si), (17) ;=I 2 The minimizer function U reads [4, I]: u ( t )  = t - $'(t) .  (18) The minimizer function x relevant to J( . ,  s) reads X ( S )  = H-I (2.4Ty + p G T s ) ,  where H = 2AT.4+ 4 GTG.  (19) 3. EXPERIMENTAL RESULTS In this section, we first consider restoring a two-dimensional image (128 x 128) and use this example to compare the per- formance of the additive and the multiplicative forms of the HQ regularization. In Figures 1 and 2, we display the orig- inal signal and the observed blurred and noisy version. We consider a spatial-invariant blurring process, and therefore the corresponding blurring matrix A is a Toeplitz-like ma- trix [14]. Since A is Toeplitz-like and G is the discretization matrix of the first-order differentiation operator, the coeffi- cient matrix ATA+BGTG in the additive form can he diag- onalized by the cosine transform matrix. It follows that the computational complexity required for solving (19) at each HQ iteration is O(n' log n) operations for an n-by-n image. However, for the multiplicative form, the coefficient matrix is ATA + pGTdiag(s)G and it cannot he diagonalized by the fast transform matrix even A and G have Toeplitz struc- tures. We employ conjugate gradient methods (inner itera- tions) to solve such linear systems. The computational com- plexity required at each inner iterations is O(n ' )  operations. The restored images using a cost-function of the form (I)-@) where $ is a Huber function (3) are displayed in Fig- ures 3 and 4 for the additive and the multiplicative forms respectively (a  = P = 1). The stopping criterion of the HQ iterations is l l f ( k )  - f ( k - 1 ) l 1 2 / / 1 g / 1 2  < 1 x 10V3. All the computations are done using MATLAB. Visually, two re- stored images using the additive and the multiplicative form are almost the same. In Table 1, we compare the perfor- mance of the HQ iterations using the additive and the multi- plicative forms for different parameters cy and @. We see that VI - 586 Fin. 1. The original image. Fig. 3. The restored image using the additive form. (relative omnr = n 1 7 A n j  Fig. 2. The blurred and noisy image. (relative error = 14.5) the multiplicative form of the HQ regularization is more ef- fective than the additive form in terms of the objective func- tion values and the relative errors. However, the differences are not significant. On the other hand, the computational time required by the additive form of the HQ regularization is significantly less than that by the multiplicative form. As we have men- tioned that fast cosine transform can be used to solve the linear system in the additive form, but inner iterations are re- quired to solve the linear system in the multiplicative form. Thus the additive form is more efficient. Next we consider the problem of increasing the spatial resolution of three-dimensional fMRI images [IZ]. A slice of a real image is given in Figures 5-8 using the additive and the multiplicative forms of HQ regularization (a = 0.01 and = 1). Both forms restore the images quite well. Again the fast cosine transformcan be applied to solve the correspond- ing linear system in the additive form. The inner iterations are required in the multiplicative form. Our numerical re- sults show that the additive form takes 103 seconds for the restoration, but the multiplicative form takes almost an hour for the restoration. 4. CONCLUDING REMARKS We performed a numerical comparison of the two forms of HQ regularization, multiplicative and additive. The ob- tained results clearly stipulate that the additive form is more Fig. 4. The restored image using the multivlicative form (relative error = 0.1239) a = 1 and0 = 1 Number of iterations required Objective function value Relative error of the restored image CPU time required (seconds) a = 0.5 and 0 = 1 Number of iterations required Objective function value Relative error of the restored image CPU time reauired lseconds) a = 0.25 and 0 = 1 Number of iterations reauired Objective function value Relative error of the restored image CPU time required (seconds) a = l a n d p = 0 . 5  Number of iterations required Objective function value Relative error of the restored imaee ICPU time required (seconds) a = l a n d p = 2  Number of iterations required Objective function value Relative error of the restored image CPU time reauired lseconds) additive I multiplicative 4 I 4 259.3545 1 259.1436 0.12408 0.12399 156.1 158.0 additive multiplicative '"+_- 159.3879 I 145.2088 0.13388 0.12751 431.9121 468.3296 0,12243 0.12215 2.8 154.2 Table 1. Comparisons between the additive and the multi. plicative forms for Huber potential function. VI - 587 Fig. 5 .  ‘lhc origitul image. Fig. 6. The blurred and noisy image. (relative error = 5.521) attractive in terms of computational cost 5. REFERENCES [I] G. Auhen and L. Vese, “A variational method in image re- covery,” SIAM Journal on Numerical Analysis, vol. 34, no. 5 ,  pp. 1948-1979, 1997. [2] Ch. Bouman and K. Saucr. “ A  generalized Gaussian image model far edge-preserving M A P  estimation,” rEEE Trans- actions on Image Processing, vol. IP-2, no. 3, pp. 296-310, July 1993. [3] M. Black and A. Rangarajan, “On the unification of line processes, outlier rejection, and robust statistics with appli- cations to early vision,’’ Inrernarionol Journal of Compurer Vision, vol. 19, no. I ,  pp. 57-91, 1996. [4] P. Charhonnier, L. Blanc-Fkraud, G. Auhert, and M. Barlaud, “Two deterministic half-quadratic regularization algorithms for computed imaging,” in Proceedings oJlEEE K I P ,  1994, vol. 2 ,  pp. 168-172. [5] P. Charhonnier, L. Blanc-FAraud, G. Auhert, and M. Barlaud, “Deterministic edge-preserving regularization in computed imaging,” IEEE Transactions on Image Processing, vol. 6, no. 2, pp. 298-31 I ,  Feb. 1997. [6]  L. Cohen, “Auxiliary variables and 2-step iterative algo- rithms in computer vision problems,” Journal of Marhemar- ical Imaging and Vision, vol. 6, no. 3, pp. 376-386, Mar. 1996. [7] D. Geman and G. Reynolds, “Constrained restoration and recovery of discontinuities,“ IEEE Transactions on Parrern Fig. 7. The restored image using the additive form. (relative rrrnr = 0.04X4X) Fig. 8. The restored image using the multiplicative form. (relative error = 0.04756) Analysis and Machine Intelligence, vol. PAMI-14, no. 3, pp. 367-383, Mar. 1992. [XI D .Geman and C. Yang, “Nonlinear image recovery with half-quadratic regularization,” IEEE Transactions on Image Processing, vol. IP-4. no. 7, pp. 932-946, July 1995. [9] R. Tyrell Rockafellar, Convex Analysis, Princeton University [IO] I. Idier, “Convex half-quadratic criteria and auxiliary inter- acting variables for image restoration,” IEEE Transactions on Image Processing, vol. 10, no. 7, pp. 1001-1009, July 2001. 11 I] P. Komprohst, R. Deriche, and G. Auhert, “Image sequence analysis via partial differential equations,’’ Journal oJMarh- emarical Imaging and Vision, vol. 11. no. I ,  pp. 5-26, Oct. 1999. [12] P. Komprohst, et al. “Superresolution in MRI and its influ- ence in statistical analuysis,” INRIA, July 2002. [I31 M. Nikolova and M. Ng, ‘Comparison of the main forms of half-quadratic regularization” Proceedings of IEEE ICIP, 2002. Press. 1970. [I41 M. Ng, R. Chan, and W. Tang. A fast algorithm for dehlur- ring models with neumann boundary conditions. SIAM Jour- nal on ScienriJic Computing, 2135 1-866, 2000. VI - 588 

Half-quadratic regularization for MRI image restoration

http://hub.hku.hk/bitstream/10722/46609/1/76681.pdf

Half-quadratic regularization for MRI image restoration

Abstract

Similar works

Full text

Available Versions

HKU Scholars Hub