A Modified ADMM Approach for Blind Deblurring of Color Medical Images

A Modified ADMM Approach for Blind Deblurring of Color Medical Images

1K Reddy Lahari

1(M.Tech), PG Scholar,

Sreenivasa Institute Of Technology &Management

Studies, Chittoor.

2V Vijaya kishore

2M.Tech,(P.hd), Professor,

Sreenivasa Institute Of Technology &Management

Studies, Chittoor.

Abstract-An efficient and flexible tool that optimizes inverse problems like image reconstruction and image restoration is alternating direction method of multipliers (ADMM)with the knowledge of blur,later ADMM is modified to perform blind image deblurring (BID)of unknown blur on original image by some function of regularization. But in real world deblurring problems the prior knowledge of blurring filter has significant importance. In this paper estimates of image and blurring operator are obtained by considering significant image edges. An ADMM iteration criterion is based on whiteness measurement which includes auto variance, auto correlation and auto covariance estimation. In the proposed approach RGB medical image of heart under different degradation conditions are considered to estimate the performance and to bring a conclusions on the degradation process and restoration on composite and component processing of the input image.

Keywords: Blind imagedeblurring,alternatingdirection method of multipliers (ADMM), non-cyclic deconvolution, Whiteness measurement, Frames-based analysis, synthesis.

1. INTRODUCTION

   In general, extract approximated original image at receiver is image restoration. In restoration main concept is deconvolving. Deconvolution is usedwhen we know an observed image blurred by a known blur kernel and degraded by an additive Gaussian noise. We use a matrix notation so that the convolution of image x with kernel R is written as Rx, where R is a block Toeplitz(diagonal-constant)matrix with Toeplitz blocksand x is taken as a column vector got by stacking all columns of the image to one long vector. In this notation, our observation model can be written as y=Rx+n, where n is a Gaussian noise of having variance 2 andmatrix R has more columns than rows because observation y includes only pixels not influenced by the unknown area outside of image x.

   Deconvolutionis usually viewed from the probabilistic viewpoint as a maximum a posterioriprobability problem, i.e., we look for image x (Original image) with the highest posterior probability, given an estimate of image prior probability distribution p(x).The main problem of the implementation in the Fourier domain isthe introduction of boundary artifacts caused by the fact that R isnot circulate (periodic).

   To solve linear inverse problem in images, we use Image restoration/reconstruction, datingback to the 1960s [3]. In this class of problems, a noisy indirectobservation, of an original image, is modeled as

   = + (1)

   Where B is the matrix representation of the direct operator and noise (n).In the particular case of image deconvolutionB is the matrixrepresentation of a convolution operator. This type of model describes well several physical mechanisms, such as relative motion between the camera and the subject (motion blur), bad focusing(defocusing blur), or a number of other mechanisms.

   In more general image reconstruction problems, B represents some linear direct operator, such as tomographic projections(Radon transform), a partially observed (e.g., Fourier) transform, or the loss of part of the image pixels.

   Next alternative method to solve several imaging problems is Alternating direction method of multipliers (ADMM), originally proposed in the 1970s [24], Emerged recently as a flexible and efficient tool compare todenoising [23], Deblurring [2], inpainting [2], reconstruction [5], motion Segmentation [4], to mention only a few classical problems (for a comprehensive review, see [6],[7],[9],[10],[13],[17],[19],[21]). It uses of variable

   splitting, which allows straightforward treatment of various priors/regularizes [1], such as those based on frames or on total-variation (TV) [3], [6], as well as the seamless inclusion of several types of constraints (e.g., positivity) [23], [8]. ADMM is closely related to other techniques, namely the so- called Bregmfan and split Bregmanmethods [10], [12], [16],

   [18],[20] [22], [23] and DouglasRachford splitting [9],[14], [19], [21].Several ADMM-based algorithms for imaging inverse problems require, at each iteration, solving a linear system [2], [10],[23]. This fact is simultaneously a blessing and accurse. The system applicability image deconvolution algorithms can also be solved by using ADMM algorithms, fast Fourier transform (FFT)is used due to simplicity, for performing the inversion of system matrix , as long as the convolution is cyclic/periodic(or assumed to be so), thus diagonal in the discrete Fourier domain.

   Generally image blurring can be obtained by convolving original with impulse response,ADMM can be used to reduce blurring effect. Practically we discussed ADMM based conjugateGradient (CG) andMask Decoupling (MD)

   algorithms in Total-Variation (TV) and Frame Based analysis and synthesis (FC) by knowing the kernels. Here finally

   (1)

   = (2)

   (4)

   observing ISNR values for four algorithms [TV-CG, TV- .

   MD, FC-CG, FC-MD]in this FC-CG will provides better results than reaming three algorithms.

   Up to now ADMM works by knowing the kernel parameters. In proposed ADMM based method will works for without knowing kernel parameters, this process is known as blind image deconvolution (BID).In proposed approach by

   ( )

   Where p=1 + 2 + + , and let g: be defined

   as

   j

   g U = gj u j (5)

   using whiteness measurement the ADMM iteration will

   j

   j=1

   automatically stops if ISNR value reaches to highest compared to previous values. Here in proposed method one additional parameter (random) was introduced.

   The proposed approach method is also applicable for color images by splitting composite image into component image and applied proposed approach method individual component images after getting results again converting component image into composite image.

2. ALTERNATING DIRECTION METHOD OF MULTIPLIERS (ADMM)

   Consider an unconstrained optimization problem

   Where each u is a -dimensional sub-vector of u, that is,u = [(u 1 ), . . , (u j ) ]

   The definitions in the previous paragraph lead to an

   ADMMfor problem (3) with the following two key features.

   1. The fact that f (z) = 0 turns line three of Algorithm1 into an unconstrained quadratic problem. Letting be a p- dimensional positive block diagonal matrix

      = i , , i , . j , , j , . , J , , (6)

      The solution of this quadratic problem is given by

      arg min ( ) = ()1

      n

      min f + g (2) z J J

      zR

      ( j ( j j )1 j ( j j (7)

      wheref : Rn .R and g : Rp .R are convex functions, and p×n is a matrix. The ADMM algorithm for this

      j=1

      j=1

      problem(which can be derived from different perspectives, namely, augmentedLagrangian and DouglasRachford splitting) ispresented in Algorithm 1 and this can be applied for

      Conjugate Gradient (CG). Algorithm 1

      1. Initialization: set k=0, choose>0, U0 and d0

      2. Repeat

         with (j) = k (j) + k (j),where ( ), ( ) ( )forj = 1, . . .

         , J , are the sub-vectors of , uk, and dk,respectively, corresponding to the partition in (4), andthe second eqality results from the block structure ofmatrices G (in (4)) and (in (6)).

   2. The separable structure of g (as defined in (4)) allows decoupling the minimization in line four of Algorithm 1 into J independent minimizations, each of the form

   ( ) arg min g + j () 2 (8)

   1. z

      arg min f + 2

      +1

      zd j 2 2

      k+1

      z 2 k k 2

   arg min f + 2

   For j=1, 2 J, Where( ) = ( ) ( ). This

   k+1

   2 k+1 k 2

   +1

   5. k+1 k (zk+1 k+1)

   6. k k + 1

   ADMM can gives exact results in linear convergence which

   minimization defines to the so-called Moreau Proximity Operator of (see [15],[16] and reference there in) applied to ( ), thus

   holds >0 however the choice of strongly effects the

   +1

   j

   min 1

   2

   convergence speed and replace scalar by diagonal-matrix

   = (1, 2 p )

   arg

   x2 2

   +

   j

   (9)

   Consider two or more (J 2) functions. Eq (2) is replaced

   by

   j

   zRn

   min gj j z (3)

   j=1

   Wheregj: are proper, closed, convex functions,

   From some functions, the Moreau proximity Operators are known in closed form [15]; a well-known case is the 1 norm for which the proximity operator is the soft-threshold function.

   Soft (, ) =sign (v)max , 0 , (10)

   Algorithm 2

   1. Initialization: set k=0, choose , . . , > 0,u d0

      and( ) × are arbitrary matrices.We map this problem intoform (1) as follows: let () = 0 define matrix G as

   2. Repeat

   3. k + k

   1 j 0,

   J J At this point, one could be tempted to map (12) into (2)using

   4.

   ( ( j (j))1 ( j (j)

   the same generalized concepts.

   k+1

   j

   j=1

   j

   j=1

   5. for j=1 to J do

   6. k+1(j) proxg u ((j)k+1 (j))

   (1) m ×d , (1) = (13)

   7. k+1

   8. end

   (j)

   k

   j j

   (j) ((j)k+1 k+1 (j))

   The problem with this choice is that the matrix to be inverted in line four of Algorithm 2 would become

   9. k k + 1

   1. until stopping criterion is satisfied

      Where the sign, max, and absolute value functions are applied in component-wise fashion. The convergence of the resulting instance of ADMM (shown in Algorithm 2)

3. DECONVOLUTION WITH UNKNOWN BOUNDARIES

1. The Observation Model

   To separate observed pixels from unobserved pixels (i.e. for bounders pixels) here we are applying one spatial mask.

   + 2 1 14

   The trick used in Periodic Boundary Condition is to express this inversion in the DFT domain can nolonger be used due to presence of M, invalidating the correspondingFFT-based implementation of line four of Algorithm2. It is clear that the source of the difficulty is the productMA, which is the composition of a mask (a spatial pointwiseoperation) with a circulate matrix (a point-wise operationin the DFT domain); to sidestep this difficulty, we need todecouple these two operations, which is achieved by defining

   g : n , g = 1 2 (15)

   1

   = + (11 )

   1 2 2

   Where M {0, 1}m×n(with m < n) is a masking matrix, i.e., a matrix whose rows are a subset of the rows of an identitymatrix. Consider that A models the convolution with a blurring filter with a limited support of size (1 + 2 ) ×

   g2: d , g2 = 1 (16)

   (1) n×d , (1) = (17)

   (2) d×d , (2) = (18)

   The only change is the proximity operator of the newg1,

   proxg

   (1 + 2 ), and thatx and Ax represent square images of 1 1

   dimensions × then matrix n×n, with = ( 2 )2 , represents theremoval of a band of width l of the outermost pixels of thefull convolved image Ax.

   Problem (11) can be seen as hybrid of deconvolutionandinpainting [14], where the missing pixels constitute theunknown boundary. IfM = I, model (11) reduces to a standardperiodic deconvolution problem. Conversely, if A = I, (11)becomes a pure inpaintingproblem. Moreover, the formulation(11) can be used to model problems where not only the boundary, but also other pixels, is missing, as in standard imagein painting.

   The following subsections describe how to handle observation model (11), in the context of the ADMM-based deconvolutionalgorithms reviewed in the previous section.

   The following sub-section describes how to apply observation model (11) in the contents of the ADMM based deconvolution algorithm which describes frame based analysis and synthesis formulation.

2. Frame-Based Synthesis Formulation

   It can be describes based on two different models that are Mask Decoupling and Conjugate Gradient.

   1. Mask Decoupling (MD): Under observation model (11),the general frame-based synthesis formulation

      = argmin 2 + 2 (19)

      u 2 1 2

      = + 1 + 1 20

      Notice that, due to the special structure of M, matrix is diagonal, thus the inversion in (20) has O(n) cost, the same being true about the product , which corresponds toextending the observed image y to the size of x, by creating a boundary of zeros around it. Of course, both + 1 and can be pre-computed and then used throughout thealgorithm, as long as 1 is kept constant. We refer to this approach as mask decoupling (MD).

      In conclusion, the proposed MD-based ADMM algorithm for image deconvolution with unknown boundaries, under frame- based synthesis regularization, is Algorithm 2 with line four implemented as Periodic Boundary Condition[1] and the proximity operators inline six given by (20) . We refer to this algorithm FS-MD (frame synthesis mask decoupling). As with the periodic BC [1], the leading cost is O(n log n) per iteration. Finally, convergence of the FS-MD algorithm is guaranteed by the following proposition (the proof of which is similar to that of Proposition 2).

   2. Using the Reevesorel Technique: An alternative to the approach just presented of decoupling the convolution fromthe masking operator is to use the method of [12], [20]

   zd 2 1

   = arg

   min 1/2 2 + z 12

   totackle the inversion (14). Following [12], notice that

   = (21)

   Where B contains the rows of AW that are missing from MAW(recall that the rows of M are a subset of those of an identity matrix) and S is a permutation matrix that puts these missing rows in their original positions in AW. Noticing that

   =

   = + (22)

   (S is a permutation matrixthus = ), the inverse of (14)can be written (with = 2 1) as

   + 1 = +

   = (23)

   Where the second equality results from using the Sherman- MorrisonWoodbury matrix inversion identity, after defining = + .Since A is circulate, C can be efficiently computed via FFT, as explained[1].Theinversion Periodic Bounder Condition cannot be computed via FFT however, its dimension corresponds to the number of unknownboundary pixels (number of rows in B), usually much smallerthan image itself. As in [12], [20], we use the CG algorithm tosolve the corresponding system; we confirmed experimentally that (as in [20]) taking only one CG iteration (initialized with the estimate of the previous outer iteration) yields the fastest convergence, without degrading the final result. Thus, in our experiments, we use a single CG iteration per outer iteration. Approximately solving line four of Algorithm 2 via one(or even more) iterations of the CG algorithm, rather than anFFT-based exact solution, makes convergence more difficult to analyze, so wewill not present a formal proof. In a related problem [23], it was shown experimentally that if the iterative solvers used to implement the minimizations defining theADMM steps are warm-started (i.e., initialized with the values from the previous outer iteration), then the error sequenceskand k, for k = 0, 1, 2,are absolutely hummableas required by Theorem 1. Finally, we refer to this algorithm asFS-CG (frame synthesis conjugate gradient).

3. Frame-Based Analysis Formulation

1. Mask Decoupling (MD): The frame-based analysis formulation corresponding to observation model Periodic Boundary Condition [1] is

   refer to this algorithm asFA-MD (frame analysis mask decoupling). The computational

   cost of the algorithm is ( )per iteration, as in theperiodic BC [1] case. Convergence of FA-MD is addressed by the following proposition

   The algorithm FA-MD (i.e., Algorithm 2 with the definition in (15), and some equations same as in Periodic Boundary Condition with line four computed and proximity operators in line six as given in (20) converges to solution of (24)

2. Using the Reevesorel Technique: Consider Problem(24) and map into Periodic BC [1] case

   (1) n×n , (1) = (25)

   The matrix inverse computed in line four of Algorithm 2 isnow (with = 2 1)

   + = + (26)

   Which can no longer be computed as in Periodic Boundary Condition [1], since matrix MAis not circulating? Using the same steps as in (21)(23), with Areplacing AW and + , leads to

   + = (27)

   Since A is circulating, + can be efficientlycomputed via FFT, as Periodic BC [1]. As in the synthesis case, the inverse in (27) is computed approximately by taking one (warm-started) CG iterationwe refer to the resultingAlgorithm as FA-CG (frame analysis conjugate gradient).

   FA-CG will gives best ISNR [1] values compared to other methods i.e., TV-CG, TV-MD, FA-MD, next we are discussing FA-CG with BID and unknown boundaries.

   IV BID using ADMM

   BIDmethod is also applicable for real images. In proposed method five different blurs (out-of-focus, linear motion, uniform, Gaussian and random) effect are explained atBSNR (90dB). The reason why we concentrate on large blur is that effect of the boundary is very evident in this case.

   By minimizing the cost (Mean Square Error) unctionthe image x and the blurring operator H (equivalently, the filter

   = arg

   min 1/2 2 +

   (24)

   1. are estimated.

      zn

      2 1

      1. Image estimating by using ADMM

         Following the MD approach introduced for the synthesis formulation, we map Problem (24) into the form (2), by usingg1 as defined in (15), and keeping H(1), H(2), and g2 as in the periodic BC case (16), (17), and (18), respectively. The only difference in the resulting instance of Algorithm 2is the use of the proximity operator of the new g1, as given in(20). In conclusion, the proposed ADMM-based algorithm for imagedeconvolution with unknown boundaries, under framebasedanalysis regularization, is simply Algorithm 2, and the proximity operators in line six given by (20). We

         The unconstrained formulation of image estimate update problem of Algorithm 1 (line 3) is defend as

         1

         2 2

         = 2 2 + 2 (28)

         =1

         And in constrained formulation (3) by letting J= m+1, and

         ( ) = , = 1, . . , (29)

         ( +1) = (30)

         2

         = 2 , = 1, . , (31)

         1

         1. BID method on color images

            The proposed method also applicable for color imageand

            +1 +1 =

            +1 2 (32)

            2

            2

            obtained ISNR values by applying different types of impulse responses, and again same color image can be split in to

            The main steps of the Algorithm 2 are those in line 4 and

            6.

            Line 4 can be written as

            k+1

            m

            ( T + + + (33)

            individual RGB component images, finding ISNR values by applying same impulse responses, which were applied on the RGB image and averaging all individual component vales. Finally comparing ISNR values of actual image with average ISNR value of individual component images.

            V METHODOLOGY

            j=1

            1. In the proposed method RGB CT heart image is taken as

               To implement line 6 of algorithm we need to solve two Proximity Operators (PO)

               g1 1 2 2

               prox v = arg q + 1 v x 2

               2

               = v-shrink (v, ,q) (34)

               prox v = arg 1 2 + 1 v x 2

               input.

            2. The image is blurred with five different blurring which is obtained by convolving the blurring function and the input image.

            3. A spatial mask is applied to prevent the unobserved pixels which may be missed during the transmission

               g1 1

               2 2 2

               process in real time application.

            4. Tight frame approach is applied to analyze the image and

               = + + (35) The proximity operator in (35) can be easily computed: is the extension of y:n by zero- padding, is a binary diagonal matrix, with zeros

               corresponding to the unobserved boundary pixels.

      2. Blur Estimating by using ADMM

         The estimate of blur problem in line 4 of algorithm 1 can be written in unconstrained formula as

         min 1/2 2 + + ()

         perform the synthesis.

         1. By using Conjugate Gradient (CG) minimizing the cost function (Mean Square Error) to get best ISNR and best restored image.

         2. The iteration process is new approach in the method, which is based on the whiteness measurement. The process of restoration will stop when the image gets best ISNR value.

         3. In addition to this iteration process the filter that responsible for the degradation of the input image is also found by Curvlets and CG methods, which is the BID.

            h 2 8. The ISNR values and the number of iterations are

            and in constrained form (3), with J=2, (1) = , (2) = and

            tabulated.

            1 1 =

            1 1 2 ,

            2

            2

            2 2 = + 2 (36)

            1. The entire approach mentioned in above steps is applied to the same image by decomposing the image in to component R, G, and B images.

            2. The individual values are tabulated and the average of

               The resulting instance of Algorithm 2 involves (in line 4) the

               inversion of the matrix(1) + (2) which can be efficiently computed in the DFT domain, using the FFT. Concerning the two (J=2) Proximity Operator (PO) in line 6, the PO has exactly the same form as (35).

               Orthogonal projection on that v-shrink is

               prox() = prox(v) = (v) (37)

               the component values is obtained along with the composite image from individual processed output images.

            3. The validation of the approach is done by comparing the ISNR values of approach on original RGB and individual component images.

         VIEXPERMENTAL RESULTS

         (2) (2)

      3. WHITENESS MEASURES

+

+

The performance analysis and validation of the proposed method is performed by considering CT image of

Whiteness measures of the residual have been previously used to assess model accuracy however; those works are on very different areas. The best new stopping criterion based on whiteness measures was successfully tested on a large resolution images, leading to an average decrease of 3% of the ISNR compared to what is obtained by stopping the algorithm at the maximum ISNR (something that, of course, cannot be done in practice, as it requires the original image).

heart. Fig1 shows original RGB image of heart CT and the resultant image obtained by applying the procedure mentioned in methodology, the obtained INR values are tabulated in Table1 for different distortions. The original CT RGB image is suppurated into component images as shown in Fig 2, the corresponding resultant images obtained by applying proposed method on individual component images are also shown in Fig 2 and their corresponding ISNR values are tabulated in Table 2, Table 3, and Table 4. The resultant images of individual R, G and B components are combined into composite image shown in Fig 3. The individual ISNR

values obtained are averaged for different distortions as shown in Table 5. Table 6 shows a comparison of ISNR values obtained by applying the proposed approach on original RGB CT heart image and average ofindividually processed images of components original RGB CT heart image.

Fig 1: original (RGB) image and result image Table 1:ISNR values for RGB image

Fig 2: component images and there resultant images

Fig 3: combined image of individual component images

Table 2: ISNR values for red component image

Table 3: ISNR values for Green component image

Table 4: ISNR values for Blue component image

Table 5: average ISNR values for above individual component image

Table 6: comparison between original image ISNR values with average of the individual color component ISNR values.

VII CONCLUSION AND FUTURE WORK:

In general deconvolution of degraded image need to have prior knowledge about degradation function, but by using BID deconvolution of degraded image can be performed without having prior knowledge about the causes for degradation. By using BID estimating the both image and degradation function can be done. In BID ADMM iteration criteria is based on whiteness measurement. Here BID method is applied on original image and individual color components of the color image and finally compared ISNR values between original image and average ISNR values of the individual color components.Hence stating that transmitting a color image by split into individual component images will give better ISNR values than transmitting a color image directly

REFERENCES

1. M. Almeida and M. A. T. Figueiredo, Deconvolving Images with Unknown Boundaries Using the Alternating Direction Method of Multipliers,IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 22, NO. 8, AUGUST 2013.

2. M. V. Afonso, J. M. Bioucas-Dias, and M. A. T. Figueiredo, An augmented Lagrangian approach to linear inverse problems with compound regularization, in Proc. 17th IEEE Intern. Conf. Image Process.ICIP, Sep. 2010, pp. 41694172.

3. M. V. Afonso, J. M. Bioucas-Dias, and M. A. T. Figueiredo, Fast image recovery using variable splitting and constrained optimization, IEEE Trans. Image Process., vol. 19, no. 9, pp. 23452356, Sep. 2010.

4. M. V. Afonso, J. M. Bioucas-Dias, and M. A. T. Figueiredo, An augmented Lagrangian approach to the constrained optimization formulation of imaging inverse problems, IEEE Trans. Image Process., vol. 20, no. 3, pp. 681695, Mar. 2011.

5. M. V. Afonso, J. M. Bioucas-Dias, and M. A. T. Figueiredo, Non-cyclic deconvolution using an augmented Lagrangian method, in Proc. IEEE Int. Conf. Comput. Tool (EUROCON), Apr. 2011, pp. 14.

1. J. Bioucas-Dias and M. Figueiredo, A new TwIST: Two-step iterative shrinkage/thresholding algorithms for image restoration, IEEE Trans. Image Process., vol. 16, no. 12, pp. 29923004, Dec. 2007.

2. S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, Distributed optimization ad statistical learning via the alternating direction method of multipliers, Found. Trends Mach. Learn., vol. 3, no. 1, pp. 1122, 2011.

3. J.-F. Cai, S. Osher, and Z. Shen, Split Bregman methods and frame based image restoration, Multiscale Model. Simul., vol. 8, no. 2, pp. 337369, 2009.

4. A. Chambolle, An algorithm for total variation minimization and applications, J. Math.Imag. Vis., vol. 20, nos. 12, pp. 8997, Jan. Mar. 2004.

5. T. F. Chan, A. M. Yip, and F. E. Park, Simultaneous total variation image inpainting and blind deconvolution, Int. J. Imag. Syst. Technol., vol. 15, no. 1, pp. 92102, Jul. 2005.

6. P. L. Combettes and J.-C.Pesquet, A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery, IEEE J. Sel. Topics Signal Process., vol. 1, no. 4, pp. 564574, Dec. 2007.

7. P. L. Combettes and V. Wajs, Signal recovery by proximal forward backward splitting, SIAM J. Multiscale Model. Simul., vol. 4, no. 4, pp. 11681200, 2005.

8. I. Daubechies, M. Defrise, and C. De Mol, An iterative thresholding algorithm for linear inverse problems with a sparsity constraint, Commun. Pure Appl. Math., vol. 57, no. 11, pp. 14131457, 2004.

9. W. Deng and W. Yin, On the global and linear convergence of the generalized alternating direction method of multipliers, Dept. Computational Appl. Math., Rice Univ., Houston, TX, USA, Tech. Rep. 1214, 2012.

10. M. Donatelli, C. Estatico, A. Martinelli, and S. Serra-Capizzano, Improved image deblurring with anti-reflective boundary conditions and re-blurring, Inverse Problems, vol. 22, no. 6, pp. 20352053, Oct. 2006.

11. M. Elad, P. Milanfar, and R. Rubinstein, Analysis versus synthesis in signal priors, Inverse Problems, vol. 23, pp. 947968, Apr. 2007.

12. E. Esser, Applications of Lagrangian based alternating direction methods and connections to split Bregman, Dept. CAM, Univ. California, Los Angeles, CA, USA, Tech. Rep. 09-31, 2009.

13. M. A. T. Figueiredo and R. Nowak, An EM algorithm for wavelet based image restoration, IEEE Trans. Image Process., vol. 12, no. 8, pp. 906916, Aug. 2003.

14. M. A. T. Figueiredo and J. M. Bioucas-Dias, Restoration of poissonian images using alternating direction optimization, IEEE Trans. Image Process., vol. 19, no. 12, pp. 31333145, Dec. 2010.

15. M. A. T. Figueiredo and R. D. Nowak, An EM algorithm for waveletbased image restoration, IEEE Trans. Image rocess., vol. 12, no. 8, pp. 906916, Aug. 2003.

16. T. Goldstein and S. Osher, The split Bregman method for l1- regularized problems, SIAM J. Imag.Sci., vol. 2, pp. 323343, 2009.

17. R. Liu and J. Jia, Reducing boundary artifacts in image deconvolution, in Proc. 15th IEEE Int. Conf. Image Process., Oct. 2008, pp. 505508.

18. S. Mallat, A Wavelet Tour of Signal Processing. San Diego, CA, USA: Academic Press, 2009