From Bayesian point of view, we get the joint regularization for the estimation of image and PSF. The resulting method attempts to minimize double cost functions subject to constraints such as non-negativity conditions of the image and energy preservation of PSFs. The objective of the convergence is to minimize double cost functions by combing the energy function of the estimation of PSFs and images. Following a Bayesian paradigm, the ideal image f, the PSF h and an observed imageg fulfill
P(f, h|g) = p(g|f, h)P(f, h)
p(g) ∝p(g|f, h)P(f, h) (5.41)
Based on this form, our goal is to find the optimal ˆf and ˆhthat maximize the posteriorp(f, h|g).
J(f|h, g) =−log{p(g|f, h)P(f)}andJ(h|f, g) =−log{p(g|f, h)P(h)}express that the energy cost J is equivalent to the negative log-likelihood of the data.
The proposed double variational regularization functional in a Bayesian framework in the BV spaces is formulated according to
J( ˆf ,ˆh) = Z
Ω
(g−ˆh∗fˆ)2dxdy
| {z }
fidelity Term +λ
Z
Ω
φ(x, Dfˆ)dxdy
| {z }
image Smoothing +β
Z
Ω
|∇ˆh|dxdy
| {z }
psf Smoothing
(5.42)
wheredxdy=dxdy. The estimates of the ideal imagef and the PSFhare denoted by ˆf and ˆh, respectively, which can be iteratively alternating minimized (AM) [289]. The image smoothing term is a variable exponent, nonlinear diffusion term [48]. The PSF smoothing term represents the regularization of blur kernels.
5.3.1 Alternating Minimization of PSF and Image Energy
During the numerical computation, we compute ∇f instead of Df. Furthermore, we remove the singularity when |∇f|= 0, by approximatingJ(f) by Jε(f) withε >0 a small parameter.
Although the most common algorithm has been based on the lagged-diffusivity technique [39], [43], [250] using an iterative procedure, we can also use it for solving the denoising and deblurring respectively. Therefore, the data-driven variant exponent, linear growth div operator becomes,
div
φ
x, q
ε2+|∇fˆ|2
= div(φ(x,∇fˆ)) (5.43)
As we have discussed in the previous chapter, the scale problem between the minimization of the PSF and the image is avoided using the alternate minimization approach. We propose to solve the joint regularization equations in an alternate minimization approach with decreased complexity. The formulation is derived from Eq. (5.42) in the following:
f∈BVinf(Ω)
Jε( ˆf ,ˆh) = 12 Z
Ω
(g−ˆh∗fˆ)2dxdy+λ Z
Ω
φε(x,∇f)dxdyˆ +β Z
Ω
(∇h)dxdyˆ (5.44)
5 Data-Driven Regularization for Variational Image Restoration in the BV Space
This equation is in strictly convexity. We can also solve this equation in the alternating mini-mization (AM) algorithm for theaugmentedenergy using two partial differential equations with respect to the imagef and P SF. Furthermore, the lagged-diffusivity algorithm corresponds to exactly the AM algorithm for the augmented energy withf(n)→H(n)→f(n+1). This algorithm is used the continuous section for achieving blur identification, and data-driven image restora-tion. The two equations derived from Eq. (5.44) are using finite differences which approximate the flow of the Euler-Lagrange equation associated with it,
∂Jε/∂fˆ = α1ˆh(−x,−y)∗(ˆh∗fˆ−g)−λdiv
φ(x,∇fˆ)
(5.45)
∂Jε/∂hˆ = α2fˆ(−x,−y)∗( ˆf ∗ˆh−g)−β∇ˆh·div ∇ˆh
|∇ˆh|
!
(5.46)
In the alternate minimization, blur identification including deconvolution, and image restoration including denoising are processed alternatingly for the estimation of the image and the PSF.
The partially recovered PSF is the prior for the next iterative image restoration and vice versa.
The algorithm is described in the following:
Initialization:g(x) =g(x), h0(x) is random numbers while(nmse > threshold
(1). nth it.fˆn(x) = arg min( ˆfn|ˆhn−1, g),fix ˆhn−1(x)
(2). (n+ 1)th it.hˆn+1 = arg min(ˆhn+1|fˆn, g), fixfˆn(x), h(x)>0 end
The data-driven diffusion term in Eq.5.46 is numerically approximated in the following, div(φ(x,∇fˆ)) =
|∇f|ˆp(x)−2
| {z }
Coefficient
[(p(x)−1)∆ ˆf
| {z }
IsotropicTerm
+ (2−p(x))|∇fˆ|div( ∇fˆ
|∇f|ˆ)
| {z }
CurvatureTerm
+∇p· ∇fˆlog|∇fˆ|
| {z }
HyperbolicTerm
] with
p(x) =
( q(x)≡1 +1+k|∇G1
σ∗I(x)|2, |∇fˆ|< β
1, |∇fˆ| ≥β
(5.47)
We indicate with div the divergence operator, and with∇and ∆ respectively the gradient and Laplacian operators, with respect to the space variables. The Neumann boundary condition [2] ∂N∂fˆ(x, t) = 0 on ∂Ω×[0, T] and the initial condition ˆf(x,0) = f0(x) = g in Ω are used, wherenis the direction perpendicular to the boundary, g is the observed image. The numerical implementation of the nonlinear diffusion operator is based oncentral differencesfor coefficient and the isotropic term, minmod scheme for the curvature term, and upwind finite difference scheme in the seminal work of Osher and Sethian for curve evolution [213] of the hyperbolic term based on the hyperbolic conservation laws. We use here the minmod function, in order to reduce the oscillations and to get the correct values of derivatives in the case of local maxima and minima.
The image is restored by denoising in the process of edge-driven image diffusion as well as de-blurring in the process of image deconvolution. Firstly, the chosen variable exponent of p(x)
(a) (b) Figure 5.9: Strength ofp(x) in the Lena image. (a) Strength ofp(x) between [1,2] in the Lena image.
(b) Strength ofp(x) is shown in a cropped image with size [50,50].
is based on the computation of gradient edges in the image. In homogeneous flat regions, the differences of intensity between neighboring pixels are small; then the gradient ∇Gσ become smaller (p(x) → 2). The isotropic diffusion operator (Laplace) is used in such regions. In non-homogeneous regions (near edge or discontinuity), the anisotropic diffusion filter is chosen continuously based on the gradient values (1 < p(x) < 2) of edges. The reason is that the discrete chosen anisotropic operators will hamper the recovery of edges [177]. Secondly, the nonlinear diffusion operator for piecewise image smoothing is processed during image deconvo-lution based on a previously estimated PSF. Finally, coupling estimation of PSF (deconvodeconvo-lution) and estimation of image (edge-driven piecewise smoothing) are alternately optimized applying a stopping criteria. Hence, over-regularization or under-regularization is avoided by pixels at the boundary of the restored image.
5.3.2 Self-Adjusting Regularization Parameter
We have classified the regularization parametersλin three different levels. Here, we present the method for the selection of window-based regularization parameters λw (window w based λw, 1st level). When the window size is amplified to the size of the input image, λbecomes a scale regularization parameter for the whole image (2nd level). If we fixλfor the whole process, then the selection of regularization parameter is conducted on the level of one fixedλ for the whole process (3nd level). We assume that the noise is approximated by an additive white Gaussian noise with standard deviation σ to construct a window-based local variance estimation. Then we focus on the adjustment of parameter λand the operators in the smoothing term φ. These two computed components can be prior knowledge for preserving discontinuities and detailed textures during the image restoration. The Eq. 5.44 can be formulated in the following,
arg min Z
Ω
φ(x, Df(x, y))dxdy subject to Z
Ω
(g−hf)2dxdy
where the noise is Gaussian distributed with varianceσ2. λcan be a Lagrange multiplier in the following form,
λ= 1 σ2|Ω|
Z
Ω
div [φ(x, Df(x, y))] (g−hf)dxdy (5.48)
5 Data-Driven Regularization for Variational Image Restoration in the BV Space
Figure 5.10: a|b. (a) Computedλ∈[0.012,0.028] values in sampling windows for the image with size [160, 160]. (b) Zoom in (a) for showing the distribution of the regularization parametersλw.
λ is a regularization parameter controlling the “balance” between the fidelity term and the penalty term. The underlying assumption of this functional satisfies kfkBV(Ω) = kfkL1(Ω)+ T V(f) in theBV space. The distributed derivative |Df|generates an approximation of input
“cartoon model”and oscillation model [164]. Therefore, this process preserves discontinuities during the elimination of oscillatory noise. We note that the term R
Ω(g−hf)2dxdy is the power of the residue. Therefore, there exists a relationship among the non-oscillatory sketch
“cartoon model” [173], [33], oscillation model [164] and the reduced power of the original image with some proportional measure.
We formulate the local varianceLw(i, j) in a given windoww based on an input image.
Lw(i, j) = 1
|Ω|
Z
Ω
[fw(i, j)−E(fw)]2w(i, j)didj (5.49)
where w(i, j) is a normalized and symmetric small window, E(fw) is the expected value with respect to the windoww(i, j) on the size of the estimated image f in each iteration. The local variance in a small window satisfies var(fw) = Lw(i, j). Thereby, we can write λ for a small windoww according to Euler-Lagrange equation for the variation with respect to f Therefore, the regularization equation with respect to the windows becomes
J ε(f) =X
λwLw(i, j) +Sp(f) (5.50)
where λwis a λ in a small window w. C is a constant is a small window. gw and fw is the observed image and the estimated image in a small windoww. Thus, we can easily get many λw for moving windows which can be adjusted by local variances, shown in Fig. 5.10. These λw are directly used as regularization parameters for adjusting the balance during the energy optimization. They also adjust the strength of diffusion operators for keeping more fidelity during the diffusion process. The related regularization parameters β and γ incorporate λ, while the parameterλof the fidelity term needs to be defined.
During image restoration, the parameter λ can be switched among three different levels. The window-based parameter λw and the scale-based (entire image) parameter can be adjusted to find the optimal results. Simultaneously,λthus controls the image fidelity and diffusion strength of each selected operator in an optimal manner.