AN IMPROVED PDE-CONSTRAINED OPTIMIZATION FLUID REGISTRATION FOR IMAGE MULTI-FRAME SUPER RESOLUTION

. The main idea of multi-frame super resolution (SR) algorithms is to recover a single high-resolution image from a sequence of low resolution ones of the same object. The success of the SR approaches is often related to a well registration and restoration steps. Therefore, we propose a new approach based on a partial differential equation (PDE)-constrained optimization fluid image registration and we use a fourth order PDE to treat both the registration and restoration steps that guarantee the success of SR algorithms. Since the registration step is usually a variational ill-posed model, a mathematical study is needed to check the existence of the solution to the regularized problem. Thus, we prove the existence and of the well posed fluid image registration and assure also the existence of the used second order PDE in the restoration step. The results show that the proposed method is competitive with the existing methods.


Problem formulation
In the presence of many degradation factors, the acquired images are always in a non desired resolution. In fact, the obtained images are decimated, noised and also blurred. We consider that the LR images are taken in the same conditions with one sensor. The link between the HR image (described by a vector of size [ 2 2 × 1] where is the enhancement factor) and the associated low resolution frames (represented by a vector of size [ 2 × 1]), is given by = + ∀ = 1, 2, . . . , , where : the convolution operator describing the blur of size [ 2 2 × 2 2 ].
: the decimation operator which is related to the desired resolution of the HR image of size [ 2 × 2 2 ]. : are the motions or the warp matrix of size [ 2 2 × 2 2 ], representing random transformation between the LR frames.
: is the additive Gaussian noise vector of size [ 2 × 1]. The aim of multi-frame super resolution is the recovery of an ideal HR image . Due to the complexity of the problem we split it in two procedure [15]: 1. Approximating the transformation matrix between each couple of LR images and merge the HR image with noise and blur. 2. Computing the HR image through the blurring and noisy one .
We start by the approximation of the warp matrix and also by ensuring the existence of the solution in a suitable functional space [1,9,12,27].

The construction of the warp matrix
To construct the operators , we use a fluid registration algorithm. This approach is based on the resolution of a variational minimization problem using fluid regularization term and an 2 based data term to model the optical flow constraint. For that, we used bicubic upsampled̂︀ of the LR sequence . We then selected the reference image between the upsampled imageŝ︀ that is used for the super-resolution process aŝ︀ (1 ≤ ≤ ). Through the imagê︀ , we compute the optical flow to all other input images (for ̸ = ). We then get flow vector fields : Ω −→ R 2 . To simplify the notations we set ( ) the pixel value corresponding to the th frame such as the coordinate = ( 1 , 2 ) ∈ Ω ⊂ R 2 (Ω is the domain containing all the pixels). The image fluid registration problem is formulated as ( ) = ( ( )) for = 2, . . . , and ∀ ∈ Ω. (3.1) Our goal is to find the velocity deformation that links each image to the reference one, where = + ∇ . , (3.2) is the partial time derivative operator and are the deformations between each frame. Unfortunately, this problem is ill-posed. We have therefore to choose an appropriate regularization operator . Since we know the success of the fluid registration to handle different problems in image registration [36], we propose to use it in a well-posed functional framework. The image fluid registration problem is defined as where ( , ) = ( , , ) + ( 1 ( ) + 2 ( )), such that is a solution of where represents the set of admissible transformations, while is the set of admissible velocity deformations, and is the regularisation parameter.
To solve the minimisation problem (3.3), we use the BFGS algorithm [42]. For that, we have to obtain some optimality conditions on the functional ℒ.

Necessary Optimality conditions
In this section, we formally give the optimality conditions of the optimization problem, since we need them to compute the BFGS Algorithm. Before showing the optimality conditions, we give briefly in the following the gradient of 2 .
Proposition 4.1. Let ℎ 1 be a perturbation of such that + ℎ 1 in , then we have the following results Proof. In the proof, we assume that ( ) = (∇ , + ). = 0 Ω. Then the derivative of 2 with respect of in direction ℎ 1 is given by: since ( ) = 0 and (∇ , + ). = 0 in Ω, we have: Now we give the optimality conditions of the problem (3.3) in the following theorem.
Let be a solution of the problem (3.3) with its associated state and a Lagrange multiplier, which are satisfying the following conditions: Proof. In order to compute the optimality conditions, we must define first the Lagrangian function, which is given as follows: where is a Lagrange multiplier. We know from the theorem 3.1 that the problem (3.3) has at least one solution, this means that ℒ admits , and as saddle points, which give us the following iteration: where ℎ (for = 1, 2, 3) are the derivative direction for each argument of ℒ. After developing derivation for each argument of the Lagrangian function and using the proposition above, we find and We can deduce the following

Algorithm details
We now give the details of each step of the BFGS Algorithm 1 described previously. We will now explain in detail the implementation of the operators , and . We begin by the downsampling operator which transforms the ([ 2 2 × 1]) HR image to the ([ 2 × 1]) LR sequence, where is the decimation factor. For example, if = 2, the decimation matrix is given as follows The BFGS algorithm to compute , and , respectively.

Evaluate
Set +1 := + 1 1 and +1 := + 2 2 and +1 := Determine step length +1 such as: +1 := end if 11. end while where ⊗ represents the Kronecker product, and the matrix 1 represents the 1D low pass filtering given by To construct the operators , we used the proposed fluid registration implementation (see Algorithm 1. For that, we used bicubic upsamplingŝ︀ of the LR sequence . We then selected the reference image between the upsamplings imageŝ︀ that is used for the super-resolution process aŝ︀ (1 ≤ ≤ ). Through the imagê︀ , we compute the optical flow to all other input images (for ̸ = ). We get then flow vector fields : Ω −→ R 2 . Since in the studied super-resolution model (2.1), we formulate the optical flow as a linear operator : However, taking into account the large size of this operator , its storage is difficult and also its computation is not feasible. Therefore, we use directly the fields to warp the HR image with respect to the input LR images using bicubic interpolation. Finally, to compute the blurring operator , we use a simple Gaussian kernel. The blur operator is then calculated using the kernel size 3 , where is the standard deviation. We have now the necessary ingredients to implement the Algorithm 1. For the fusion step, we use the algorithm in [25,48] to compute the blurred and noisy HR image = .
After that, we carry out the last step of the SR algorithm, which is the restoration step.

Restoration step
Since the problem of restoration is ill-posed, we have to be careful in the choice of a suitable approach for denoising and deconvolution step. The main purpose of this stage is to preserve image features and avoid the blocky effect while reducing noise and blur. Since first order regularizers have always suffered from the blocky effect [10], numerous approaches have been proposed to tackle this effect, namely, higher-order regularizations [8,31,33,34]. Another successful approach in the super resolution framework was also proposed in [22] to preserve the details of the obtained image. A more robust higher-order total variation model is proposed as an efficient solution to the blocky effect, called total generalised variation (TGV) [44]. Even if this method avoids the staircasing effect, the computational time to reach the solution is very significant. An alternative way is to use a combined approach of the first and second order regularizer proposed by Papafitsoros et al. [37] which eliminates noise and blur, avoiding the block artifacts, in less time than the TGV. In the same esprit, we introduce a firstorder operator to preserve edges and a second-order one to remove noise in smooth areas. Let us describe briefly the algorithm used in the proposed SR approach. In this paper, we use the gradient descent PDE with a Neumann boundary condition on Ω, associated to the variational model proposed in [37]. This PDE with the initial condition (0, ) = 0 (where 0 is the obtained HR image calculated by the interpolation of the LR image 1 ) is given by where the divergence operator div : (R Also, the two divergence operator div 2 : (R , with the adjointness property, is defined We give briefly the discrete version of the operators ∇ and div given by and where (6.8) Let us define the second order discrete differential operators noted ∇ 2 as and (6.12) In addition, for = ( 1 , 2 , 3 , 4 ) ∈ (︀ R )︀ 4 , we define the discrete div 2 operator as where ∇ = ∇ , ∇ = ∇ , (6.14) and To give a comprehensive form of this problem, we set the particular case where 1 ( ) = 2 ( ) = . As a result, the algorithm related to solve numerically the proposed PDE is finally given in Algorithm 1.

Experiments
In this part, our aim is to test the ability of the elaborated algorithm in the SR context. Many simulated and also real results were used to test the performance of the proposed SR method. The first part is dedicated to the evaluation of the registration part while the second and third parts concern the main proposed SR approach.

The effectiveness of the registration part
In the first experience, we construct 32 synthetic LR images from the original image of MRI Brain 1 such that: each frame is deformed by random vector fields, blurred by a Gaussian low-pass filter with size 4 × 4 and a standard deviation of 2. Then the blurred frames are down-sampled in the two directions by a factor of = 3 and a Gaussian noise was added with a standard deviation = 20. In this test, we fix the deconvolution and denoising part where we use the proposed combined first and second order regularizer. Then, we compare our registration algorithm with other competitive registration methods in the SR context, such as: SR with probabilistic optical flow (POF) [16], SR with hyper-elastic (SRHE) [29], SR with consistent flow (SRCF) [51] and also SR with diffusion registration (SRDR) proposed in [25]. The obtained SR results are shown in Figure 1 to see the efficiency of the proposed registration part. We can deduce that the proposed registration method  [16]. (f) SRCF [51]. (g) SRHE [29]. (h) SRDR [25]. (i) Our method.
gives a slightly better result compared to the other methods. For the second experiment, we consider the MRI Brain 1 following the same previous steps with more complicated random deformations between the LR images. We increase also the decimation factor = 4 and also the standard deviation = 35. The obtained HR image is depicted in Figure 2, where comparison to other SR methods is done. Clearly the obtained HR image outperforms the other ones and the registration process is better enhanced.
In the following tests, the first six experiments were simulated ones with known HR images, while the next three experiments were real data experiments. In all simulated experiments, we measure the effectiveness of our method by considering some comparisons with popular SR methods, such as bicubic interpolation, [15], SR method [49] and also the combined first order and ( + ) regularisation [24] using the same data. In the following experiments, the motion model is assumed to be the global translational model for the other methods while we use the fluid registration for our approach. In this paper, in the simulated experiments, the peak-signal-to-noise ratio (PSNR) [21] and also the structural similarity (SSIM) [45] are computed to check the quality of the recovered HR image.

Simulated Experiments
In the simulated experiments, we choose six examples in Figure 3, with different size and texture, taken from a tested benchmark images to illustrate the performance of the proposed SR method. To construct the simulated images, we follow the degradation model (2.1). The HR image was first shifted with sub-pixel displacements to produce images, then, the sequence was convoluted with a PSF and finally, zero-mean Gaussian noise was added to each frame of the sequence. For example, in the first experiment, we take the image of Boat 3a as an original image of size 512 × 512, we construct then = 50 input low-resolution frames by shifting the Figure 3a in vertical and horizontal directions, sub-sampling with a decimation factor = 4 and blurred using a Gaussian density with kernel size 3 × 2. Finally, we add a white Gaussian noise in each frame with a standard deviation = 10. We use the same thing for the Gear image, while we increase the blur kernel rate and noise for the other tests. Indeed, we use a 5 × 5 Gaussian blur kernel and a white Gaussian noise with = 30 to construct the sequence for the last four images. To increase the ability of the proposed equation to better detect edges, we choose the so-called hypersurface minimal function [3] defined for both the functions 1 and 2 by The input HR image 0 is built by a bicubic interpolation of the LR reference image 1 . Concerning the choice of the parameters, the scalar weight is chosen according to the better PSNR value for the proposed and also for the other SR methods. For instance, we choose = 0.6 for the Boat image. Concerning the convergence of the proposed algorithm, we end the execution at the first iteration with respect to the error ‖̂︀ +1 −̂︀ ‖ 1 ‖̂︀ ‖ 1 < 10 −5 . We set also = 10 −3 . reg. [15]. (d) method [49]. (e) + reg. [24]. (f) The proposed SR. reg. [15]. (d) method [49]. (e) + reg. [24]. (f) The proposed SR. reg. [15]. (d) method [49]. (e) + reg. [24]. (f) The proposed SR. reg. [15]. (d) method [49]. (e) + reg. [24]. (f) The proposed SR. reg. [15]. (d) method [49]. (e) + reg. [24]. (f) The proposed SR.
The obtained results using the proposed SR are illustrated and compared with other approaches in the Figures 4-9 for the tested image in Figure 3 respectively. Visually, we can see that the proposed model gives better restoration than the others. The difference becomes more obvious when the noise and blur increase which is the case in the Figures 6-9. The effectiveness of the proposed approach model can also be shown in the color images in Figures 8 and 9, where the obtained HR image is more clean with fewer registration errors. However, to approve the success of the proposed algorithm against noise and misregistration errors reducing, we use the PSNR criterion in the Table 1 and the SSIM measure in the Table 2 for three noise values. Knowing that the best score is in bold number, we can clearly show the efficiency of our algorithm. Also, visually, we can detect the performance of the proposed method in removing misregistration errors compared with the other methods in the smooth area, however, in the edge area there is no distinct improvement compared to + and methods.

Real Experiments
In the real experiments, three real data sequences are used to approve the proposed algorithm, are presented. We select the first ten frames in the three real data sets. The registration approach presented in [44] is used as the registration estimation method for the other methods while we use the proposed fluid registration to estimate the motion for our method. The reconstruction results of these sequences are, respectively, shown in  reg. [15]. (d) method [49]. (e) + reg. [24]. (f) The proposed SR. + reg. [24]. (c) method [49]. (d) Our method.  Figure 11. Results on the Text sequence. (a) First LR frame. (b) + reg. [24]. (c) method [49]. (d) Our method.  Figure 12. Results on the Emily sequence. (a) First LR frame. (b) + reg. [24]. (c) method [49]. (d) Our method.

Conclusion
In this paper, a novel approach of the super resolution image reconstruction problem is introduced. We presented a bilevel fluid image registration and proved the existence and uniqueness of the solution using functional analysis theory. In addition, to avoid the undesirable staircasing effect during the registration, denoising and deconvolution steps, a fourth-order PDE is proposed. To show the robustness of this approach, a set of benchmark image have been performed, and the investigated SR approach has proven its success visually and also quantitatively using two known metrics. A remaining question is about the treatment of other types of noise and blur, such as Salt& paper and Poisson noise. Another interesting point is about the degrees of the efficiency of the proposed approach with respect to the nature of the transformations between the LR frames.