## Abstract

Image denoising is important for high-quality imaging in adaptive optics. Richardson-Lucy deconvolution with total variation(TV) regularization is commonly used in image denoising. The selection of TV regularization parameter is an essential issue, yet no systematic approach has been proposed. A construction model for TV regularization parameter is proposed in this paper. It consists of four fundamental elements, the properties of which are analyzed in details. The proposed model bears generality, making it apply to different image recovery scenarios. It can achieve effective spatially adaptive image recovery, which is reflected in both noise suppression and edge preservation. Simulations are provided as validation of recovery and demonstration of convergence speed and relative mean-square error.

© 2014 Optical Society of America

## 1. Introduction

Adaptive optics is a promising technique applied to compensate wavefront distortion and de-noise acquired images. Besides high definition imaging in telescope and remote sensing equipment, it also has widespread applications, such as biomedical science, astronautics and laser atmospheric transmission. Adaptive optics basically includes research work on optical distortion compensation and image denoising. As a result of optical system design, closed loop servo bandwidth and capability of computer processing, state-of-the-art optical distortion compensation suffers from insufficient correction. Image denoising is necessary to achieve noise suppression and high definition imaging. In fields of biomedical science, astronautics and remote sensing, there is usually demanding requirement for noise suppression and definition. In the transmission channel of these fields, the essential part is the convolution of point spread function and original image. Consequently, many deconvolution algorithms are provided as approaches, such as blind deconvolution [1] (specially, iterative blind deconvolution (IBD) in [2]), Richardson-Lucy(RL) deconvolution [3, 4], etc. In practice, RL deconvolution has proved to be robust and reliable and is considered to better than IBD in deconvolution. As a result, it has widespread applications [5], and many approaches around RL deconvolution have been proposed [6, 7].

However, RL deconvolution itself bears an inherited defect due to the ill-posed nature of matrix inversion [8]. As a result, some regularized versions of RL deconvolution are brought forward as solutions to the problem. Mathematically, there exist some choices of regularization models, e.g. ridge regression, Lasso and basic pursuit(BP) denoising. In [9], Rudin *et al.* propose the total variation(TV) regularization. Because of the differential properties of gradient function, TV regularization exhibits outstanding edge preservation results [10, 11], which is of great importance in noise suppression and image enhancement. Its idea is to achieve solution with the most sparse gradient, and in theory, E. Candes *et al.* has proved in [12] that this approach is conditionally capable of recovering information exactly from incomplete observations. It allows realization of approximate inverse computation [13]. Besides image processing, TV regularization has also seen widespread applications in X-ray imaging, super-resolution, coherent imaging, computed tomography and holography [14–23].

In TV regularization, the essential issue is the regularization parameter *λ*. Much of previous work has adopted TV regularization for image denoising and recovery, and it is related to the TV regularization parameter [24, 25]. However, no systematic approach of its selection has been proposed up to now. Therefore, we propose a systematic model to fill the gap in this paper. The scope of our work lies in RL deconvolution, and we propose a construction model for TV regularization parameter. This model consists of four fundamental elements, i.e. constructor function, edge indicator, intensity factor and shaping factor. The definitions, properties and effects of the four fundamental elements are analyzed in details in the following sections. The TV regularization parameter selection is generalized into the combination of the four fundamental elements, and this model is adequate for a wide variety of applications. Simulations of image recovery by means of the proposed model are illustrated to validate the construction model and demonstrate the convergence speed and relative mean-square error.

## 2. Mathematical description of the model

#### 2.1. Introduction of RL deconvolution and TV regularization

In adaptive optics, the isoplanatic linear imaging system is expressed as

where**,**

*i***,**

*h***and**

*o***are the observed image, the point spread function, the original image and noise, respectively. The symbol * refers to convolution. In Eq. (1),**

*n***∈**

*x***denotes pixels**

*X***in the image space**

*x***, generally represented as (**

*X**x*,

*y*) in a two-dimensional Cartesian coordinate system. Moreover, the point spread function

**has a normalized imposition, i.e. ∑**

*h*

_{x}**(**

*h***) = 1.**

*x*In the imaging system described in Eq. (1), consider the Poissonian noise of ** n**. Thus, the conditional probability is represented as follows:

In RL deconvolution, Bayesian theory is adopted to solve the maximum a posteriori(MAP) and minimize the functional

In TV regularization, a regularization parameter *λ* is added and the functional ** J**(

**) becomes**

*o*When we minimize the above Eq. (4) with respect to ** o** and

**and search for zero points respectively, then apply D’Alembert convergence condition $\frac{\Vert {\mathit{o}}^{n+1}\Vert}{\Vert {\mathit{o}}^{(n)}\Vert}=1$, we obtain**

*h**n*+ 1) and (

*n*) are numbers of iterations, ‖·‖ represents

*ℓ*

_{2}-norm, div(·) refers to divergence,

*o*^{*}and

*h*^{*}are conjugates of

**and**

*o***, respectively, and Eq. (7) is for preservation of the normalization imposition of**

*h***.**

*h*#### 2.2. Approximation in weak or moderate noise cases

In weak or moderate noise cases, the intensity of ** n** is negligible or relatively small, i.e. Eq. (1) is approximate as

As a result, with conditions of Eq. (8) and the normalized imposition of ** h**, Eq. (5) is rewritten as

The regularization parameter *λ*, which is the only artificially introduced parameter, is of essence in TV regularization. In flat areas, noise contamination is the major problem, and the expression in Eq. (10) is small. Therefore, *λ* is required to be large to achieve noise suppression. On the other hand, in higher gradient areas, there exists detailed information we want to preserve, and the divergence in Eq. (10) is large. As a result, *λ* has to be small to preserve the details in higher gradient area, and the TV regularization algorithm basically degenerates to the unregularized RL deconvolution. In this way, *λ* should be a bounded function with respect to ** x**, and the ideal selection of

*λ*(

**) will achieve spatially adaptive regularization:**

*x**λ*(

**) is large in flat areas and small in higher gradient areas, and we define it as the**

*x**spatially adaptive property*of

*λ*(

**).**

*x*#### 2.3. Proposition of our model and properties of λ(**x**) and fundamental elements

According to the discussion above, we propose the construction model for TV regularization parameter *λ*(** x**) as follows:

*f*is

*constructor function*,

*EI*(

**) is**

*x**edge indicator*,

*λ*

_{0}is

*intensity factor*and

*β*is

*shaping f actor*. The two factors of

*λ*

_{0}and

*β*are positive constants, and they get their names from their vertical and horizontal effects on the graph of

*λ*(

**), respectively.**

*x**f*,

*EI*(

**),**

*x**λ*

_{0}and

*β*are

*fundamental elements*of the regularization parameter construction model, and the combination of the four fundamental elements makes the construction process systematic and convenient.

Next, we will analyze properties of *λ*(** x**) and the four fundamental elements. As discussed before,

*λ*(

**) is a bounded function, therefore, let**

*x**m*be the minimum (or infimum) and let

*M*be the maximum (or supremum). In this way,

*λ*(

**) ∈ (**

*x**m*,

*M*]. The edge indicator

*EI*(

**) :**

*x***→ ℝ**

*X*^{+}∪ {0} is a non-negative function defined in the image space

**. It is used to distinguish the image “edgy degree”. Its value increases with the increase of image “edgy degree”, i.e.**

*X**EI*(

**) is small in flat areas and large in higher gradient areas. In order to adapt all edge indicators**

*x**EI*(

**) satisfying the aforementioned conditions to the proposed model, the domain of**

*x**f*is set as [0, +∞). If we consider the spatially adaptive property of

*λ*(

**) and the property of**

*x**EI*(

**) comprehensively, it is obvious that the constructor function**

*x**f*is a monotonic decreasing function. Its domain is [0, +∞) and its range is (

*m*,

*M*]. Furthermore,

*f*(

*x*) should be convex when

*x*is sufficiently large, since it is a monotonic decreasing function with the domain of [0, +∞) and the bounded range of (

*m*,

*M*].

The above Eq. (9) also indicates information of convergence, which provides instruction for the choices of *m* and *M*. It is clearly seen that if we want Eq. (9) to converge faster, *λ*(** x**) should be small, i.e. the minimum

*m*has to be near zero, which is related to detail preservation in higher gradient areas. For the value of the maximum

*M*, we can set it empirically, as long as it is sufficiently large for noise suppression in flat areas. The speed of

*f*(

*x*) → 0 should also be fast when

*x*→ +∞, which translates as

*λ*(

**) fast arrival of zero point for higher gradient areas. This guarantees fast convergence of Eq. (9).**

*x*Then, we provide some qualitative discussion of the intensity factor *λ*_{0} and the shaping factor *β*. Since the regularization parameter *λ*(** x**) being too small will lead to insufficient noise suppression and it being too large will cause blur of image detailed information,

*λ*

_{0}being too small and/or

*β*being too large will fail to suffice noise suppression, and

*λ*

_{0}being too large and/or

*β*being too small will blur the image details. As their names suggest,

*λ*

_{0}changes the overall intensity of

*λ*(

**) and**

*x**β*shapes

*λ*(

**) horizontally. Therefore, the values of**

*x**λ*

_{0}and

*β*should be chosen comprehensively.

In brief, the essential properties of four fundamental elements of the construction model proposed in this paper include:

- Constructor function
*f*(Property 1): The constructor function*f*is a monotonic decreasing function, and it should be convex when its argument is sufficiently large. Its domain is [0, +∞) and it ranges from (*m*,*M*], where*m*is near zero and*M*is an empirical value large enough to suppress noise in flat areas. - Edge indicator
*EI*() (Property 2): The edge indicator*x**EI*() is a non-negative function defined in the image space*x*. It increases with the increase of image “edgy degree”, i.e. it is small in flat areas and large in higher gradient areas.*X* - Intensity factor
*λ*_{0}and shaping factor*β*: The parameters of*λ*_{0}and*β*are positive scalar constants.

## 3. Simulation

In order to verify the image recovery effects of our proposed model, we will construct TV regularization parameter *λ*(** x**) with the four fundamental elements as in Eq. (11), and simulation results of image recovery with the regularization parameter are provided.

To begin with, we preset the bounds of *m* and *M* as *m* = 0 and *M* = 1, and we choose some constructor functions according to Property 1. We choose *f* in the category of elementary function, and some typical choices are *f*_{1}(*x*) = 1/(1 + *x*), *f*_{2}(*x*) = 2/*π* · arc cot(x), *f*_{3}(*x*) = 1/(1 + *x*^{2}) and *f*_{4}(*x*) = exp(−*x*). The typicality of *f*_{1} ∼ *f*_{4} lies in the fact that there are six types of fundamental elementary functions, i.e. constants, exponentials, logarithms, *n*-th roots, trigonometric functions and inverse trigonometric functions, and the choices of *f*_{1} ∼ *f*_{4} cover all possible types. Graphs of *f*_{1} ∼ *f*_{4} are shown in Fig. 1.

Besides, we want to demonstrate the performance of TV regularization with a fixed *λ* for comparison with our proposed model, therefore we choose *f*_{0}(*x*) = 1, since *f*_{0} indicates a fixed value of *λ*(** x**) =

*λ*

_{0}. The images in the following Fig. 2 and Fig. 4 are all greyscale images in the range of 0 to 255.

#### 3.1. Recovery of a real noisy image

We will recover a real moon image blurred by noise, the size of which is 218 × 195. For the edge indicator *EI*(** x**), we choose a very intuitive one, the

*ℓ*

_{2}-norm of gradient (first-order derivative) of the image, i.e. $\mathit{EI}(\mathit{x})=\Vert \nabla \mathit{o}\Vert (\mathit{x})=\sqrt{{\left(\frac{\partial \mathit{o}}{\partial x}\right)}^{2}+{\left(\frac{\partial \mathit{o}}{\partial y}\right)}^{2}}(x,y)$. It obviously satisfies Property 2. The two parameters of

*λ*

_{0}and

*β*are constants set as

*β*= 0.05 and

*λ*

_{0}= 0.02, respectively. The recovery process goes through 400 iterations. In Fig. 2, the blurred moon image and its five recovered versions corresponding to constructor functions

*f*

_{0}∼

*f*

_{4}are demonstrated. Figure 2(a) is the blurred moon image and Figs. 2(b)–2(f) are images recovered by constructor functions

*f*

_{0}∼

*f*

_{4}, respectively. It is clearly seen that Fig. 2(b) suffers from serious staircase effect, which is typical in fixed TV regularization parameter RL deconvolution. Moreover, the recovery is more and more effective from Fig. 2(c) to Fig. 2(f) and we can see more delicate details in the image.

To demonstrate the convergence of iterations, we provide the convergence results of Figs. 2(b)–2(e) in Fig. 3. Due to the D’Alembert convergence condition
$\frac{\Vert {\mathit{o}}^{(n+1)}\Vert}{\Vert {\mathit{o}}^{(n)}\Vert}=1$ we apply in Eq. (5), we use the expression of relative error of successive iterations
$\frac{\Vert {\mathit{o}}^{(n+1)}-{\mathit{o}}^{(n)}\Vert}{\Vert {\mathit{o}}^{(n)}\Vert}$ to examine convergence of iterations. As it is seen, *f*_{0} does not converge at all. All constructor functions of *f*_{1} ∼ *f*_{4} realize fast convergence, and the iterations converge faster from *f*_{1} to *f*_{4}: the relative error of successive iterations becomes smaller for the same number of iterations. The convergence results of *f*_{4} stand out since the error is the smallest among the four constructor functions and it generally decreases when the number of iterations increases.

#### 3.2. Recovery of a simulated noisy image

An original circuit image with the size of 256 × 256, its blurred version and recovered versions are shown in Fig. 4 to further demonstrate recovery results of simulated noisy image. We choose to adopt a more sophisticated edge indicator, which is a modified *D*(** x**) from [26], as

*EI*(

**). The modification lies in the Hessian matrix of the image, where we add a Gaussian filter**

*x**G*and makes convolution of

_{σ}*G*and each element of the Hessian matrix. (The size of Gaussian filter is 3 × 3 and the deviation

_{σ}*σ*is 0.5.) The functions of

*f*

_{0}∼

*f*

_{4}shown in Fig. 1 are still employed as constructor functions. For the intensity factor and the shaping factor, we maintain

*λ*

_{0}= 0.0105 for all four constructor functions, and set

*β*= 2.0 for

*f*

_{1}∼

*f*

_{3}and

*β*= 0.85 for

*f*

_{4}. The fourth

*β*is smaller because ${f}_{4}(x)=\text{exp}(-x)=1/{\sum}_{n=0}^{+\infty}{x}^{n}/n!$, and in Fig. 1 it is clearly seen

*f*

_{4}decreases faster with

*x*increasing for the same

*β*as

*f*

_{1}∼

*f*

_{3}. The recovery process goes through 300 iterations. We can see Figs. 4(d)–4(g) exhibit visually the same recovery, which are all better than Fig. 4(c) since Fig. 4(c) suffers from insufficient luminance and obvious noise blur.

Figure 5 provides the convergence results and relative mean-square error. In Fig. 5(a), the same expression of relative error of successive iterations as in Fig. 3 is employed to verify convergence. To examine the enhancement of image recovery, we adopt the relative mean-square error of the recovered image with respect to the original image
$\frac{{\Vert \widehat{\mathit{o}}-\mathit{o}\Vert}^{2}}{{\Vert \mathit{o}\Vert}^{2}}$, where ** o** and

**respectively represent the original image and the recovered image. Relative mean-square error results are shown in Fig. 5(b). It is obvious that**

*ô**f*

_{0}fails to converge and the relative mean-square error of

*f*

_{0}is not well suppressed. As we can see, the relative mean-square error results of

*f*

_{1}∼

*f*

_{4}are almost the same, while the convergence results of

*f*

_{1}∼

*f*

_{4}differ. If we take both convergence and relative mean-square error into account,

*f*

_{2}and

*f*

_{3}perform better, which show both high converge speed and low relative mean-square error.

It is worth mentioning that our construction model is generalized, and this generality is guaranteed by the aforementioned analyses in Section 2 (esp. Section 2.3). Videlicet, the constructor function *f* (*x*) and the edge indicator *EI*(** x**) satisfying Property 1 and Property 2 (including but not limited to examples provided in the simulations) with appropriate values of

*λ*

_{0}and

*β*will achieve proper recovery results.

## 4. Conclusion

In conclusion, this paper presents a construction model for TV regularization parameter, which is expressed as *λ*(** x**) =

*λ*

_{0}·

*f*(

*β*·

*EI*(

**)). It consists of the constructor function**

*x**f*, the edge indicator

*EI*(

**), the intensity factor**

*x**λ*

_{0}and the shaping factor

*β*as the four fundamental elements. The properties of the fundamental elements are provided in Property 1 and Property 2. Image recovery by means of the construction model becomes spatially adaptive, which is reflected as noise suppression in flat areas and detail preservation in higher gradient areas. Simulation curves of convergence and relative mean-square error illustrate high convergence speed and low recovery error of this model. Moreover, the proposed model bears generality, which offers a systematic way to construct TV regularization parameter, making it applicable in biomedical science, astronautics, remote sensing and other applications.

## Acknowledgments

This work was supported by the National Basic Research Program of China (Grant No. 2013CB329201) and the National High Technology Research and Development Program of China (Grant No. 2013AA013601). The authors are more than grateful to the handling editor and the reviewer for their valuable suggestions to this paper.

## References and links

**1. **A. V. Oppenheim, R. W. Schafer, and T. G. Stockham Jr., “Nonlinear filtering of multiplied and convolved signals,” Audio and Electroacoustics, IEEE Transactions on , **16**(3), 437–466 (1968). [CrossRef]

**2. **G. R. Ayers and J. C. Dainty, “Iterative blind deconvolution method and its applications,” Opt. Lett. **13**(7), 547–549 (1988). [CrossRef] [PubMed]

**3. **W. H. Richardson, “Bayesian-based iterative method of image restoration,” J. Opt. Soc. Am. **62**(1), 55–59 (1972). [CrossRef]

**4. **L. B. Lucy, “An iterative technique for the rectification of observed distributions,” Astron. J. **79**, 745 (1974). [CrossRef]

**5. **G. Krishnamurthi, C. Y. Wang, G. Steyer, and D. L. Wilson, “Removal of subsurface fluorescence in cryo-imaging using deconvolution,” Opt. Expr. **18**(21), 22,324–22,338 (2010). [CrossRef]

**6. **D. A. Fish, A. M. Brinicombe, E. R. Pike, and J. G. Walker, “Blind deconvolution by means of the Richardson-Lucy algorithm,” J. Opt. Soc. Am. A **12**(1), 58–65 (1995). [CrossRef]

**7. **D. S. C. Biggs and M. Andrews, “Acceleration of iterative image restoration algorithms,” Appl. Opt. **36**(8), 1766–1775 (1997). [CrossRef] [PubMed]

**8. **D. A. Hope and S. M. Jefferies, “Compact multiframe blind deconvolution,” Opt. Lett. **36**(6), 867–869 (2011). [CrossRef] [PubMed]

**9. **L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D **60**(1), 259–268 (1992). [CrossRef]

**10. **E. Vera, P. Meza, and S. Torres, “Total variation approach for adaptive nonuniformity correction in focal-plane arrays,” Opt. Lett. **36**(2), 172–174 (2011). [CrossRef] [PubMed]

**11. **M. Freiberger, C. Clason, and H. Scharfetter, “Total variation regularization for nonlinear fluorescence tomography with an augmented Lagrangian splitting approach,” Appl. Opt. **49**(19), 3741–3747 (2010). [CrossRef] [PubMed]

**12. **E. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inform. Theory **52**(2), 489–509 (2006). [CrossRef]

**13. **A. Kostenko, K. J. Batenburg, H. Suhonen, S. E. Offerman, and L. J. van Vliet, “Phase retrieval in in-line x-ray phase-contrast imaging based on total variation minimization,” Opt. Expr. **21**(1), 710–723 (2013). [CrossRef]

**14. **E. Y. Sidky, M. A. Anastasio, and X. Pan, “Image reconstruction exploiting object sparsity in boundary-enhanced x-ray phase-contrast tomography,” Opt. Expr. **18**(10), 10,404–10,422 (2010). [CrossRef]

**15. **A. Kostenko, K. J. Batenburg, A. King, S. E. Offerman, and L. J. van Vliet, “Total variation minimization approach in in-line x-ray phase-contrast tomography,” Opt. Expr. **21**(10), 12,185–12,196 (2013). [CrossRef]

**16. **J. I. Sperl, D. Bequé, G. P. Kudielka, K. Mahdi, P. M. Edic, and C. Cozzini, “A Fourier-domain algorithm for total-variation regularized phase retrieval in differential X-ray phase contrast imaging,” Opt. Expr. **22**(1), 450–462 (2014). [CrossRef]

**17. **A. Chambolle, “An algorithm for total variation minimization and applications,” J. Math. Imaging Vision **20**(1–2), 89–97 (2004). [CrossRef]

**18. **J. Dahl, P. C. Hansen, S. H. Jensen, and T. L. Jensen, “Algorithms and software for total variation image reconstruction via first-order methods,” Num. Alg. **53**(1), 67–92 (2010). [CrossRef]

**19. **M. Cetin, W. Karl, and A. Willsky, “Edge-preserving image reconstruction for coherent imaging applications,” in *Proceedings of ICIP 2002 International Conference on Image Processing* (Rochester, 2002). [CrossRef]

**20. **D. Brady, K. Choi, D. Marks, R. Horisaki, and S. Lim, “Compressive holography,” Opt. Expr. **17**(15), 13,040–13,049 (2009). [CrossRef]

**21. **M. Marim, M. Atlan, E. Angelini, and J. Olivo-Marin, “Compressed sensing with off-axis frequency-shifting holography,” Opt. Lett. **35**(6), 871–873 (2010). [CrossRef] [PubMed]

**22. **A. Bourquard, N. Pavillon, E. Bostan, C. Depeursinge, and M. Unser, “A practical inverse-problem approach to digital holographic reconstruction,” Opt. Expr. **21**(3), 3417–3433 (2013). [CrossRef]

**23. **E. Y. Sidky and X. Pan, “Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization,” Phys. Med. Biol. **53**(17), 4777–4807 (2008). [CrossRef] [PubMed]

**24. **L. M. Mugnier, J.-M. Conan, T. Fusco, and V. Michau, “Joint maximum a posteriori estimation of object and PSF for turbulence-degraded images,” in *Bayesian Inference for Inverse Problems*, A. Mohammad-Djafari, ed., Proc. SPIE3459, 50–61 (1998). [CrossRef]

**25. **N. Dey, L. Blanc-Feraud, C. Zimmer, Z. Kam, P. Roux, J. C. Olivo-Marin, and J. Zerubia, “Richardson-Lucy Algorithm with Total Variation Regularization for 3D Confocal Microscope Deconvolution,” Microsc. Res. Tech. **69**(4), 260–266 (2006). [CrossRef] [PubMed]

**26. **H. Tian, H. Cai, J. Lai, and X. Xu, “Effective image noise removal based on difference eigenvalue,” in *Proceedings of ICIP 2011 International Conference on Image Processing* (Brussels, 2011). [CrossRef]