Mixed Higher Order Variational Model for Image Recovery

A novel mixed higher order regularizer involving the first and second degree image derivatives is proposed in this paper. Using spectral decomposition, we reformulate the new regularizer as a weighted L 1 -L 2 mixed norm of image derivatives. Due to the equivalent formulation of the proposed regularizer, an efficient fast projected gradient algorithm combined with monotone fast iterative shrinkage thresholding, called, FPG-MFISTA, is designed to solve the resulting variational image recovery problems under majorization-minimization framework. Finally, we demonstrate the effectiveness of the proposed regularization scheme by the experimental comparisonswith total variation (TV) scheme, nonlocal TV scheme, and current second degreemethods. Specifically, the proposed approach achieves better results than related state-of-the-art methods in terms of peak signal to ratio (PSNR) and restoration quality.


Introduction
During image acquisition, the observed images are mainly degraded by blurring and random noise, which are tightly linked to imaging system.In order to attenuate the undesired effects, image recovery plays an integral part of modern imaging sciences in such areas, including medical imaging, microscopy, astronomical imaging, remote sensing, and photography [1][2][3][4][5].
Image recovery which amounts to estimating a desired image from its degraded measurements is mathematically an inverse problem [6].To cope with the ill posedness of image recovery, the standard approach is to use image a priori to constrain the solutions so that many variational regularization approaches have been successfully proposed for image recovery.As we know, one of the most popular regularizers is the total variation (TV) seminorm which was firstly introduced by Rudin et al. for image denoising [7] and turned out to be very efficient for preserving image edges.Thus, TV has been widely and effectively used in many other image applications [8][9][10][11] and provides state-ofthe-art result.However, TV-based reconstruction methods often give rise to the undesired staircase effect, since TV favors piecewise constant solutions.Thus, the main challenge for image recovery is to reduce the staircase effect while preserving edges sharpness as far as possible.As is known to us, piecewise-vanishing second degree derivatives yield piecewise-linear solutions which better fit smooth intensity changes.To alleviate the staircase effect, several regularizers involving second degree image derivatives were introduced mainly for image denoising over past few years [12][13][14][15] and provide better performance than TV.Recently, many novel higher order regularization methods which constitute valid extensions of TV are effectively introduced for more image processing tasks [16][17][18] so that they minimize the staircase effect and provide state-of-the-art results.However, the second-order regularizers [16][17][18] which are the mixed norms of second degree directional image derivatives tend to suppress larger second degree directional derivatives thus leading to blurring across the sharp edges.Thus, the reconstructed images with better preserved edges and less staircase effect are more desirable.To this end, as reported in [19][20][21][22][23], the regularizers that combine first and second degree image derivatives have been proposed so as to provide better edge preservation while reducing the staircase effect in smooth regions of image; thus, we are motivated to investigate the regularizers involving derivatives of different degrees in this paper.
In this paper, our main contributions are twofolds.First, we reinterpret the isotropic second degree total variation (ISDTV) regularizer [16] which is the  1 - 2 mixed norms of the second degree directional derivatives and give its equivalent formulation.We then propose a novel mixed higher order regularizer involving the first and second degree derivatives of the image and reformulate the new regularizer as a weighted  1 - 2 mixed norm of image derivatives in spectral decomposition framework.The new regularizer which inherits the desirable properties of TV turns out to be well suited for image recovery.Second, for the proposed variational image recovery model, a fast projected gradient algorithm combined with monotone fast iterative shrinkage thresholding (FPG-MFISTA) which is derived from majorizationminimization (MM) framework has been utilized to solve the minimization problem under additional convex constraints.Finally, we compare our proposed schemes with TV [24], nonlocal TV [25], and current second degree methods [18,26], which currently provide state-of-the-art results.
The rest of the paper is organized as follows.In Section 2, we first present the mathematic formulation of image recovery problem and review the ISDTV regularizer.We then propose a novel regularizer at the end of Section 2. In Section 3, we propose our variational regularization-based image recovery model and give a detailed description of the corresponding minimization algorithm.In Section 4, we assess the performance of our approach on both natural and biomedical images and provide systematic comparisons with TV, nonlocal TV, and current second degree methods.Finally, we conclude the paper in Section 5.

Mixed Higher Order Total Variation for Image Recovery
2.1.Regularization for Image Recovery Problem.We first consider the recovery of a continuously differentiable image  : Ω → R from its degraded observation , where Ω ⊂ R 2 is the spatial support of the image.The most commonly used image observation model can be linearly formulated as where  is a linear blur operator and  is Gaussian white noise with standard deviation .As we know, the recovery of  from its observation  is an ill-posed problem, due to the fact that the blur operator  is ill-conditioned.To stabilize the recovery process, some available a priori knowledge about the desired image is utilized.Thus, a common variational regularization approach is employed to formulate image recovery as an optimization problem: where the first term is known as data fidelity term, which measures the consistency between estimation and observation, while the second one is called regularization term which provides some a priori knowledge of the desired image, and the regularization parameter  ≥ 0 provides a balance between the two terms.

Isotropic Second Degree Total Variation and Equivalent
Formulation.Hu and Jacob have introduced the isotropic higher degree total variation (I-HDTV) scheme [16], which inherits the desirable properties of TV, such as, convexity and invariance to rotations and translations.As detailed in [16], the first degree directional derivative of  along the unit vector s 0 () = (cos , sin ) T at coordinate (, ) denoted by  ,1 (, ) is defined as where g 0 (, ) denotes the gradient of  at coordinate (, ) and (⋅) T denotes transpose operator.Specially, as proved in [16], the standard TV is reinterpreted as the  1 - 2 mixed norms of the first degree directional derivative, formulated as where ‖ ⋅ ‖ 2 denotes the Euclidean norm.Now, we consider the isotropic second degree total variation (ISDTV) of image , which is defined as where  ,2 (, ) is the second degree directional derivative of image  at coordinate (, ), defined as with Based on the work of [16], we derive an equivalent formulation of the integrand of (5) in following proposition, whose proof is provided in Appendix A.
Based on the definitions in (10), we define the new regularizer as Definition (11) leads to a convex regularizer which is also homogeneous, rotation, and translation invariant.Moreover, we have namely, the regularizer in (11) involves the first and second degree directional derivatives of the image which are related to TV and ISDTV, respectively, and takes full advantage of the ability to preserve edges and minimize the staircase effect.Thus, the resulting regularizer constitutes a valid extension of TV and ISDTV.For simplicity, we term the new regularizer in (11) as mixed higher order total variation (MHOTV).
Based on the following proposition whose proof is given in Appendix B, we give an equivalent formulation of the MHOTV regularizer.(11) is equal to ‖RG(, )‖ 2 , where ‖ ⋅ ‖ 2 stands for the Euclidean norm.
Therefore, according to the differential operator V defined in Table 1, we further express our MHOTV regularizer as With the equivalent formulations of the ISDTV and MHOTV regularizers defined in (9) and (13), which can be also understood as weighted  1 - 2 mixed norm of image derivatives, thus the two equivalent formulations are more preferable for the minimization algorithm we will propose in Section 3.

FPG-MFISTA: The Efficient Algorithm for the Proposed Model
Using the ISDTV and MHOTV regularizers described in ( 9) and ( 13), we present a unified variational regularization image recovery model to reconstruct the desired image from its noisy observation, which is formulated as where the last expression in (14) means that the regularization term can be either ISDTV or MHOTV regularizer.

Discrete Formulation for ISDTV-and MHOTV-Based
Image Recovery Model.Now we consider the discrete formulation of proposed model (14).In order to simplify our analysis, we assume that the image with size  ×  is stacked in a vector of size  =  × .Then, the objective function in ( 14) can be discretized as where A ∈ R × is the blurring matrix and g, f ∈ R  are the vectorized versions of the observed image and the image to be recovered, respectively.To discretize the differential operators defined in Table 1, we assume reflexive boundary conditions for the images and use forward finite differences to approximate the partial derivatives [27].
To simplify our analysis, we will introduce some auxiliary definitions and notations as follows.
First, we define operators U : R  → R 3× and V : R  → R 5× to refer to the discrete versions of the operators  and  defined in Table 1, respectively.For simplicity, we denote space   = R  × = {u = (u 1 , u 2 , . . ., u  ) | u  ∈ R  , ∀ = 1, 2, . . ., } and equip the space   with the inner product ⟨⋅, ⋅⟩   and norm ‖ ⋅ ‖   .Let u = (u 1 , u 2 , . . ., u  ) ∈   and k = (k 1 , k 2 , . . ., k  ) ∈   ; we define them as Next, let us define the adjoint operators of U and V which are the discrete operators U T :  3 → R  and V T :  5 → R  , respectively.Thus, let p ∈  3 , q ∈  5 and f ∈ R  ; we have where ⟨⋅, ⋅⟩ 2 stands for the inner product in Euclidean space R  .Furthermore, with the  1 - 2 mixed norm of a vector field u = (u 1 , u 2 , . . ., u  ) ∈   defined as ‖u‖ 1,2 = ∑  =1 ‖u  ‖ 2 , we thus write the discrete versions of the ISDTV and MHOTV regularizers defined in ( 9) and ( 13) as where (⋅)  denotes the th column element of the corresponding vector filed.Therefore, we can write the discrete formulation of proposed model (14) as with the last expression H ∈ {WU, RV } in (20) meaning that the regularization term can be either ISDTV or MHOTV regularizer.

Majorization-Minimization (MM).
In the following, we focus on the minimization of model (20).Due to the fact that the ISDTV and MHOTV regularizers are both nondifferentiable, we solve model ( 20) by a majorizationminimization (MM) approach [28][29][30].As described in the MM framework, an iterative algorithm for solving the proposed model (20) at th iterate f () involves the following two steps: (1) constructing the surrogate functional (f; f () ) that satisfies the following properties: (2) solving the next iterate f (+1) = argmin f (f; f () ).
Under the MM framework, we find the optimal solution of model ( 20) by iteratively minimizing a sequence of surrogate functionals {(f; f () )} that upper bound the objective function (f).
Therefore, to obtain the surrogate functionals {(f; f () )}, we first upper bound the data fidelity term of (f) at th iterate f () using the following majorizer: where z = f () +(1/)A T (g−Af () ) and  is a constant which is independent of f.To make sure that (f; f () ) is valid, we choose  > ‖A T A‖ to make I−A T A positive definite.Thus, the whole resulting majorizer of (f) at th iterate f () can be written as Then, the minimization problem ( 20) is converted to iteratively minimizing (23) with regard to f, which can be formulated as As we can see, the resulting minimization problem ( 24) is much simpler and its solution can be alternatively interpreted as the solution of a denoising problem with z being the noisy observation.In the following, we will present a primal-dual approach based on the monotone fast iterative shrinkage thresholding algorithm (MFISTA) [24] to seek the optimal solution of (24).

FPG-MFISTA.
Considering that f is a digital image whose intensities are finite and bounded, we can use a convex closed set to constrain the range of image intensities.Based on that, the resulting denoising subproblem (24) can be equivalently formulated as where H ∈ {WU, RV} means that the regularization term can be either ISDTV or MHOTV regularizer,  = {f ∈ R  | f  ∈ [, ], ∀ = 1, . . ., } is a convex closed set which enforces bounded constraints on the solution, and   (f) is the indicator function of  that takes the value 0 for f ∈  and ∞ otherwise.
Next, we present a primal-dual approach which can efficiently compute the solution of the constrained problem (25).Without loss of generality, here we just describe all the steps of the minimization algorithm for model (25) involving the proposed MHOTV regularizer, which corresponds to the case that H = RV.Therefore, by denoting  = /, the corresponding MHOTV based denoising subproblem can be further rewritten as As referred in [31], ‖ ⋅ ‖ ∞,2 is the dual norm of ‖ ⋅ ‖ 1,2 , and then the  1 - 2 mixed norm of the vector field u = (u 1 , u 2 , . . ., u  ) ∈   can be rewritten as where  = { = ( thus, the resulting formulation (28) makes us to solve the following minimax problem: Since the cost function (f, ) is strictly convex in f and concave in , we can exchange the order of the minimum and maximum, and the saddle-point (f (+1) ,  () ) of ( 29) is attained as reported in [32].Therefore, As we know, (f, ) can be equivalently written as Based on (30) and (31), we can obtain the saddlepoint (f (+1) ,  () ) of ( 29) by solving where   denotes the orthogonal projection onto the convex set , while  () is the maximizer of the dual problem, formulated as Since the dual function ℎ() in ( 33) is differentiable and has Lipschitz continuous gradient, as referred in [24], then we can compute the gradient of ℎ() as To accelerate the convergence, we use Nesterov's iterative method [33] to solve model (33), which is a gradient-based approach that exhibits convergence rates of one order higher than the standard gradient ascent method.Furthermore, we employ the operator   to return the orthogonal projection onto the  ∞ - 2 unit norm ball , where this operation is performed by projecting independently each component   ,  = 1, 2, . . .,  of  onto the  2 unit norm ball.
Therefore, as described above, we first find the maximizer of the dual problem (33) and then obtain the optimal solution of our primal problem ( 26) through (32).
In the same way, by using ( 28)- (33) and replacing R and V with W and U, respectively, we can obtain the solution of model (25) involving the ISDTV regularizer.
Finally, a detailed description of the fast projected gradient algorithm combined with monotone fast iterative shrinkage thresholding (FPG-MFISTA) for ISDTV and MHOTV regularization-based image recovery is given in Algorithm 1, while the corresponding denoising subproblem (24) can be efficiently solved by a fast projected gradient algorithm described in Algorithm 2.

Experimental Results
To evaluate the performance of our proposed ISDTV and MHOTV regularization schemes in the context of image deblurring, we provide systematical experimental comparisons with five state-of-the-art methods, namely, TV regularization method [24], nonlocal TV regularization method [25], Hessian Frobenius-norm and Hessian Spectral-norm regularization methods [18], and the Hessian Nuclear-norm regularization method [26].Finally, the quality of the reconstructed image is evaluated in terms of peak signal-to-noise ratio (PSNR), which is measured in decibels (dB) and defined as where f ∈ R  and f (0) ∈ R  are the vectorized versions of the reconstructed image and the original image, respectively.

Experiment Setting.
For the image deblurring experiments, our test images contain six natural images and six biomedical cell images which are shown in Figures 1 and 2, respectively, where the cell images are part of the biomedical image database [34], and they are all converted to grayscale.
In our experiments, the intensities of all the test images are initially normalized to the range [0, 1], and we use three kinds of point spread function (PSF) to produce the blurred images, which are a 9 × 9 Gaussian PSF of standard deviation   = 6, a 9 × 9 uniform one, and a 19 × 19 motion one, respectively.Furthermore, the Gaussian noise of two different levels corresponding to a blurred SNR (BSNR) of {20, 30} dB is added to generate the eventually degraded images, and the BSNR is defined as where var(Af) is the variance of blurred image and   is the standard deviation of Gaussian noise.In addition, we constrain the intensities of image recovered by all the seven methods to lie in the convex closed set  = {f ∈ R  | f  ∈ [0, 1], ∀ = 1, . . ., }.For the sake of fair comparisons, the regularization parameter  for each regularization method is chosen to give the best PSNR performance.Finally, we set  = 10 −5 ,  = 100,  0 = 10 −3 , and  = 10 in our experiments.

Recovery on Natural Images.
In Table 2, we provide systematic comparative results for all the regularization schemes on the natural images under different degradation conditions.As we can see from Table 2, regarding comparisons among the seven regularization schemes, the best PSNR results, on average, are achieved for the proposed MHOTV regularizer, the nonlocal TV regularizer provides better results than TV regularizer for most of the cases, and ISDTV regularizer outperforms TV and Hessian Spectralnorm regularizers for most of the cases, while the Hessian Frobenius-norm regularizer always performs similar results with ISDTV scheme, and the PSNR results of Hessian Nuclear-norm regularizer are almost always a little better than that of Hessian Spectral-norm and Frobenius-norm    regularizers.Thus, we can conclude that our MHOTV-based regularization scheme is proven to be a valid approach for image recovery.Furthermore, we show the representative Fingerprint, Face, and Peppers deblurring examples in Figures 3, 4, and 5, which are degraded by Gaussian blurring, uniform blurring, and motion blurring with different Gaussian noise levels, respectively.
From these figures, we can clearly see that the images recovered by TV regularization method suffer from heavy blocky artifacts appearing in the smoother regions of the image and result in some loss of details, the nonlocal TV regularization method preserves the sharp edges but smoothes some moderate gradient areas, the Hessian Spectral-norm, Hessian Frobenius-norm, Hessian Nuclearnorm, and ISDTV regularization methods are very competitive, and they can also preserve some sharp edges and minimize the staircase effect, while our MHOTV regularization method provides the best image quality which can effectively minimize the staircase effect and better preserve edges sharpness.For example, in Figure 3, since the Fingerprint image is mainly composed of smooth transitions and ridge features, in order to highlight the differences, we present the details of the corresponding zooming versions of the region (marked by the black box in Figure 3(a)) for all the seven methods in Figures 3(c)-3(i).Comparing the results of these seven methods, we observe that the TV solution shown in Figure 3(c) results in some loss of fine ridge details and produces the staircase effect, the nonlocal TV solution shown in Figure 3(d) smoothes some ridge details, and the remaining five methods all preserve fine ridge details and edges better than TV and nonlocal TV methods, while our MHOTV solution provides the best performance in terms of both PSNR and image quality.

Recovery on Biomedical Cell Images.
In this section, we mainly talk about the performance of the proposed MHOTV regularization scheme on biomedical images, which are closely related to the practical applications.To evaluate the effectiveness of our approach, we provide deblurring experiments for all the methods on the set of biomedical cell images, and their comparative results are given in Table 3.
With the PSNR comparisons in Table 3, we can clearly observe that the proposed MHOTV regularization approach provides the best PSNR results for most of the cases, the ISDTV method always performs better than nonlocal TV, TV, and Hessian Spectral-norm methods, and the nonlocal TV method provides the better results than TV method for most of the cases, while the Hessian Frobenius-norm method always provides the comparatively similar results to ISDTV method, and the PSNR results of Hessian Nuclear-norm method are almost always a little better than that of Hessian Spectral-norm and Frobenius-norm methods.Meanwhile, we show the representative Fluorescent cell 3 and CIL 12293 deblurring examples in Figures 6 and 7, respectively.While carefully inspecting the recovered images shown in Figures 6 and 7, we obviously find that the advantage of our MHOTV approach consists in better preserving the ridges and filament-like features sharpness without introducing severe block artifacts, and the nonlocal TV, ISDTV, and the three Hessian-norm regularization methods can preserve some fine ridge features and alleviate the staircase effect, as opposed to TV which introduces the heavy block artifacts which tends to produce the oversmoothed image ridge and filament-like features.Therefore, we conclude that our MHOTV regularization approach can be effectively and successfully applied in the process of biomedical imaging.

Conclusion
In this paper, we have proposed a novel MHOTV regularizer which combines the first and second degree derivatives of the image.Within spectral decomposition framework, we reformulate both the ISDTV and MHOTV regularizers as weighted  1 - 2 mixed norms of image derivatives.Furthermore, we present a unified variational model based on the ISDTV and MHOTV regularizers for image recovery.Based on the equivalent formulations of ISDTV and MHOTV regularizers, we then introduce an efficient fast projected gradient algorithm combined with monotone fast iterative shrinkage thresholding (FPG-MFISTA) to solve the corresponding optimization problems under majorizationminimization framework.Finally, the performance of the proposed MHOTV regularization scheme was assessed by comparing the experimental results on both natural images and biomedical cell images against various state-of-the-art methods.Specifically, the proposed MHOTV regularization scheme performs improved results and can effectively minimize the staircase effect that is commonly met in TV scheme while better preserving edges sharpness.

B. Proof of Proposition 2
Proof.Using the rotation steerability of directional derivatives [35,36], we can simplify the integrand of (11) as

Table 2 :
PSNR (dB) results of the seven regularization methods on natural images for three blurring kernels and two different noise levels.

Table 3 :
PSNR (dB) results of the seven regularization methods on biomedical cell images for two blurring kernels and two different noise levels.