Using Kalman filter in the frequency domain for multi-frame scalable super resolution

Accepted Manuscript

Using Kalman filter in the frequency domain for multi-frame scalable super resolution Akbar Rahimi , Payman Moallem , Kamal Shahtalebi , Mehdi Momeni PII: DOI: Reference:

S0165-1684(18)30296-2 https://doi.org/10.1016/j.sigpro.2018.09.012 SIGPRO 6922

To appear in:

Signal Processing

Received date: Accepted date:

14 August 2018 8 September 2018

Please cite this article as: Akbar Rahimi , Payman Moallem , Kamal Shahtalebi , Mehdi Momeni , Using Kalman filter in the frequency domain for multi-frame scalable super resolution, Signal Processing (2018), doi: https://doi.org/10.1016/j.sigpro.2018.09.012

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Using Kalman filter in the frequency domain for multi-frame scalable super resolution Akbar Rahimi1, Payman Moallem1*, Kamal Shahtalebi1, Mehdi Momeni2 1- Department of Electrical Engineering, Faculty of Engineering, University of Isfahan, Isfahan, Iran

*Corresponding Author

CR IP T

2- Department of Surveying Engineering, Faculty of Engineering, University of Isfahan, Isfahan, Iran

Abstract: Kalman filter (KF) as a linear estimator which is used in super-resolution (SR)

AN US

problems, suffers from high computational costs and storage requirements. To gain appreciable success in the elimination of these two challenges, this paper advances a SR framework employing KF in the frequency domain, while no resort is made to any approximations or extra assumptions in the dynamic system modeling and statistical matrices.

M

Generally, previous KF-based SR methods organized the system with huge-sized matrices in the spatial domain, following which they tried to reduce the system dimension using

ED

approximation and/or limitation on point spread function (PSF). In this study, first, several

PT

small-dimension dynamic systems are separately made in the frequency domain supporting space-invariant PSFs of an arbitrary form and size. Then, the acquired small-dimension KF

CE

estimators are applied rather than the traditional huge-dimension one. These will greatly reduce computational complexity, decrease storage requirements allowing parallel

AC

implementation as well. Furthermore, our proposed SR framework can be used to produce high resolution image of an expedient size, that is, a scalable SR. Experimental results with both simulated and real world sequences indicate that our proposed framework works more effectively than the other compared methods, especially in fine details restoration.

Keywords: Super-resolution, Kalman filter, frequency domain 1

ACCEPTED MANUSCRIPT

1-Introduction The higher the resolution of an image, the more the details that emerge in observation. So a digital image with higher resolution is always more desirable in such wide applications as remote sensing, military, medical diagnosis, traffic monitoring than traditional vision imaging [1]. However, the available array imaging sensor is not dense enough due to the physical

CR IP T

limitations and forbidding cost constraints. The multi-frame super-resolution (SR) reconstruction can greatly assist in extracting the requisite information from a set of low resolution (LR) degraded images by creating a high resolution (HR) image [2]. The LR images are captured at different perspectives of the same scene. Generally, at least one of

AN US

sub-pixel translations among the LR frames is considered a good assumption in the very area of real applications.

The multi-frame image super resolution framework was first introduced by Tsai and Huang

M

[3] in the frequency domain. They derived a linear system equation in accordance with the

ED

relation between LR images and the desired HR image. They also assumed that the frames were band-limited having sub-pixel translations among them. Then, Tekalp et al. [4]

PT

considered observation noise and spatial blurring in developing their modeling formulation. Later, Kim et al. [5] provided another extension of the frequency domain approaches based

CE

on weighted least squares formulation and Tikhonov regularization. To reduce the memory

AC

requirements and the computational costs, a discrete cosine transforms (DCT)-based method was proposed by Rhee and Kang [6]. Frequency domain algorithms typically establish use of the relationship between the desired HR image and the observed LR images based on a simple theory, requiring short computation time [1]. On the other hand, a variety of spatial domain methods have been developed incorporating non-uniform interpolation [7], iteration back-projection (IBP) [8], projection onto convex sets (POCS) [9], HR reconstruction without primary registration [10], Kalman filter-based 2

ACCEPTED MANUSCRIPT methods, and regularized methods. It should be noted that most of the recent research works in the spatial-based SR problems focus on the regularized frameworks because of their effectiveness and flexibility [1]. The regularized methods estimate the desired HR image by minimizing an objective function which includes two main terms: Data fidelity and regularization. While the first one measures the closeness between the reconstructed HR

resulting in a robust solution [11,12,13,14,15,16,17,18].

CR IP T

image and the captured LR images, the second one is employed to regularize the problem

Generally, tuning regularized sensitive parameters is a notorious challenge in regularized SR

AN US

methods. In comparison to the spatial domain methods, frequency domain approaches might have difficultly supporting complicated motion conditions and being sensitive to model errors. Nonetheless, they enjoy several major advantages including theoretical simplicity,

M

facility for parallel implementation, and needing no sensitive parameters tuning. Kalman filtering theory is widely known as an optimal estimation theory under linear

ED

constraints and Gaussian noise [19]. Some SR reconstruction methods are based on the Kalman filter theory, which methods try to estimate the HR image from a set of observed

PT

aliased noisy LR images with regard to a linear observation model. On the whole, the

CE

traditional Kalman filter-based SR method suffer from high computational complexity and huge storage requirements, which emanate from the storage of high dimension matrices and

AC

special type of operations like matrix inversion during updating procedures. Some researchers have tried to meet the major challenges inherent in Kalman filter-based SR methods. Dellaert et al. [20] proposed a simplified Kalman filter-based algorithm which neglects off-diagonal elements of covariance matrix for the computational feasibility. Farsiu and Milanfar [21] carried out a two-step SR algorithm. In the first step, they used the Kalman filter framework for fusing the data, and in the second one, they applied the de-blurring technique. Newland et al. [22] developed a modified Kalman filter using the steady state Kalman gain matrix and 3

ACCEPTED MANUSCRIPT square-form filter for imaging point spread function (PSF). Recently, Wei et al. [23] proposed point-wise Kalman filter-based method by introducing a new system of modeling and changing covariance matrix to modified block diagonal covariance matrix in order to solve the computation and storage problems. Although the developed Kalman filter-based approaches have overcome some of the existing disadvantages of the Kalman filters,

degrade the quality of the reconstructed HR image.

CR IP T

limitations on the PSF type and size, and the approximations on covariance matrix still

The present study introduces a novel Kalman filter-based SR method without any recourse to

AN US

approximations to achieve better performance. Also the method we have proposed requires no tuning parameters, supports the space-invariant PSFs of any form and size, and is capable of reconstructing the HR image with an arbitrary size and scale. Further, it reduces the computational complexity and the required storage simultaneously. As a first step, an

M

observation model is developed relating the LR images to the continuous scene. Then the

ED

discrete Fourier transform (DFT) is evolved to obtain the frequency domain observation model. Next, the low-dimension dynamic system equations are organized based on the

PT

spectrum observation relation. Finally, the Kalman filter is employed to estimate the system

CE

states of constituting the frequency components of the desired HR image. In fact, the frequency domain helps to develop separate low-dimension dynamic systems. In

AC

other words, all pixels of HR image in the spatial domain relate to all pixels of each LR image while in the spectrum domain, a few frequency components of HR image relate to each frequency component of LR image causing a separation of several dynamic systems and leading to several low-dimension Kalman filters. After this general overview, it is in order to present the format of what appears in the following parts. In Section 2, an observation model is developed in the spectrum domain. In

4

ACCEPTED MANUSCRIPT Section 3, we organize the dynamic system equations, derive corresponding Kalman filter estimators, and describe our proposed SR algorithm. The experimental results are presented in Section 4. And the final Section 5 brings this study to an end by enumerating the

CR IP T

outstanding characteristic achievements of the method we have proposed here.

2-Image observation model

An observation model which simulates the relation between the desired HR image and each of the observed LR images is an essential step in the SR image procedure. Supposing the

AN US

imaging model which universally involves warping, blurring, down sampling, and noise adding, relates each observed LR image to a desired HR image.

Similar to the most represented SR methods, only global sub-pixel translations from among

M

LR images are considered, i.e., those which adapt to the most real successive videos. Let us suppose scene which is considered as a two dimensional (2-D) continuous signal being

(

). To obtain each observed LR image for the continuous scene, as a

) is translated by

and

coming from camera shifting to a reference

PT

first step,

(

ED

represented as

point in -direction and -direction, respectively. Then the shifted continuous signal

(

CE

), assumed to be band limited, is sampled at or above Nyquist rate to define a new

version of the desired HR image. Blurring, down-sampling, and noise adding are the other

AC

steps in the observation modeling as are depicted in Fig.1. Further, we can express the image observation as: ( (

))

(1)

5

ACCEPTED MANUSCRIPT Where

are rth observed LR image and the added measurement noise,

and

respectively.

is the imaging system PSF.

and

are the sampling and down-sampling

operations, respectively. * is the 2-D discrete convolution operator. of

in -direction and -direction, respectively. Besides,

and

are also the size

is the number of available

CR IP T

frames.

rth Warped HR image Continues scene

rth Observed LR image

Optical or/and Motion blur and imaging system PSF

Down Sampling by factor d

AN US

Sampling at or above Nyquist rate

Continues translation

Noise

M

Fig.1, Observation model relating LR images to HR image in spatial domain.

ED

Eq. (1) formulates the observation model in the spatial domain for the SR problem to be expressed in the frequency domain. Accordingly, step by step, we transform the obtained

PT

spatial-based model to the spectrum domain. (

CE

As a first move, let

(FT) of continuous signal

AC

respectively while

( (

(

) be the Fourier transform

) and the corresponding shifted version represent frequency in

(

-direction and

), -direction,

. Based on the sampling theory, the DFT of the signal

)) is given by: (

Where

and √

respectively. Also

(

) and

)

(2)

;

is the DFT of the desired HR image

6

.

and

are the pixel sizes

ACCEPTED MANUSCRIPT in -direction and -direction, respectively. represent discrete frequency in

is also well-known as the SR factor(srf).

and

(

) is

-direction and

-direction, respectively.

assumed band limited. For the shifted continuous signal sample

( (

)) which is rth HR image, DFT

( (

))

↔

(

)

Where

and

, while

and

CR IP T

can be given as follow:

(3)

represent the pixel plus sub-pixel

AN US

shifting of the rth HR image in -direction and -direction, respectively. By comparing Eq.(3) relating DFT of the rth HR image to that of the

and Eq.(2), the translation function

reference HR image can be defined as follows:

(4)

M

,

ED

Now, assume that the PSF of the observing system is space invariant. Then we can write: ( (

)) ↔

is the DFT of PSF . For simplicity of notations, we set

CE

PT

Where

(5)

.

(6)

Finally, if we apply the DFT of the down sampled 2-D signal to factor , the observation

AC

equations can be completely expressed in the spectrum domain as ∑

Where

(7)

is the measurement noise at

DFT of the noise variance

∑

,

. It is easy to show that if

frequency component equaling the is a white Gaussian noise with

will also be a white Gaussian noise with variance 7

(cf.

ACCEPTED MANUSCRIPT Appendix A). Moreover, in cases where

is non-Gaussian, based on the central limit

theorem (CLT), when independent random variables are added, their proper sum tends toward a Gaussian distribution [24], which means each component of

can

be approximated by a Gaussian noise.

AN US

and which is to be exploited in our SR reconstruction method.

CR IP T

Fig. 2 illuminates our proposed observation model developed out of the previous discussion

M

Fig. 2, Observation model relating frequency components of rth LR image to that of HR

ED

image.

CE

algorithm

PT

3- Developing the dynamic systems, Kalman filters, and the proposed

It can be seen from Eq. (7) that the frequency domain observation model of the super

AC

resolution is a linear measurement. Hence, the Kalman filter estimator can be used to solve super resolution problems. To attain this objective, first the proper dynamic systems must be organized. 3-1-Dynamic systems

8

ACCEPTED MANUSCRIPT Essentially, every dynamic system can be set up through a state vector determination and establishing two equations: system and measurement equations. Based on Eq. (7), the states vector

can be defined as:

[

(

)

(

)

]

And, the measurement matrix is also defined as (

CR IP T

(8)

)

(

)

AN US

(9)

Now, based on the linear relation (7) and the defined matrix, the dynamic system equations are realized as follows:

(10)

ED

M

{

Where

represents the system modeling error assumed to be zero,

.

PT

represents the measurement error described above, with the covariance matrix for all , while

is corresponding to the FT of

AC

CE

the desired HR image at the reference position.

3-2-Kalman filter and the proposed SR algorithm Based on the developed dynamic systems through Eq.(10), Kalman filters can be derived to estimate the states being frequency components of desired HR image. There are two steps to the Kalman filter process. The first step is prediction, the second one is updating. For the

9

ACCEPTED MANUSCRIPT represented dynamic systems, these two steps can be joined and the corresponding Kalman filter should be derived as follows [19]: ( ̂

)

̂

̂

(

)

(11) (12)

and

are the Kalman gain and the covariance matrix,

respectively. Superscript

represents the matrix Hilbert transpose.

Eq. (11) provides the gain matrix at utilizing

gain

th

matrix ̂

observation and Eq. (12) updates the system state ̂ and

the

measurement

residual

AN US

Where

CR IP T

(13)

. Additionally, Eq. (13) updates the system covariance

matrix. Using these equations, the Kalman filtering process can be formulated and renewed . Afterwards, they are repeated iteratively.

M

for

Finally, the SR procedure is completed by placing each final version of ̂ to Eq.(8), and applying the inverse DFT of ̂

according

ED

to obtain the desired HR image.

PT

To constrain the required number of iterations in a way that the algorithm converges, we

‖̃

CE

consider that the latter algorithm is terminated if the following condition is satisfied: ̃

‖ ‖̃

‖

(14)

AC

Where ̃ is the estimate of the DFT of the desired HR image at the tth iteration, ‖ ‖ and denote norm L2 and the terminated threshold, respectively. Our proposed SR method can be summarized in the following algorithm scheme. It is noted that in our research the same initial value matrices

where

is considered for all covariance

is a presetting scalar.

10

ACCEPTED MANUSCRIPT

Algorithm

AN US

CR IP T

Predetermine values: noise variance, for all ̂ Initial values; , I (I is identity matrix), for all , ( ) Input frames ; Calculate DFT of the frames as 1- Repeat while terminated condition Eq.(14) not fulfilled 2-Repeat for to 3- Repeat for =0 to 4- Repeat for to Calculate Kalman gain matrix according to Eq.(11) Estimate ̂ according to Eq.(12) Update covariance matrix according to Eq.(13) End of repeat 4 End of repeat 3 End of repeat 2 Set End of repeat 1 Calculate inverse DFT of final ̂ as the desired HR image

M

4-Experimental Results

ED

In this section, we make use of both artificially created LR images and real frames to test our proposed SR algorithm, so we can evaluate its performance in comparison with other existing

PT

SR methods. Three metrics, PSNR (peak signal to noise ratio), SSIM [25] (structural similarity index measure) and FSIM [26] (feature similarity index measure) are employed to

CE

verify the quality of the reconstructed HR images in simulation cases. Then real world frames

AC

are applied to visually evaluate and compare our SR approach with other related methods presented in application. In our experiments, different PSFs selected from Table 1 are employed. Moreover, in the final experiments, the scalable ability of our proposed SR method to reconstruct HR images of different sizes and scales is represented. In order to have a fair validation of our proposed algorithm, we employ five SR algorithms, namely, two methods based on the Kalman filter, that is, video-to-video dynamic system SR method (VTV) [21] and the Kalman filter-based method (KFB) [23], and three regularized 11

ACCEPTED MANUSCRIPT SR methods which are L1 norm with BTV method (L1+BTV)[11], adaptive norms with modified BTV method (AN+MBTV) [12], and regional, spatially adaptive total variation method (RSATV) [13]. In all experiments, initial values are provided in the manner presented in Table 1. In addition,

CR IP T

regularization parameters are adjusted by several rounds of running in algorithms- VTV [21], L1+BTV [11], AN+MBTV [12], and RSATV [13]- until the best HR results are achieved. Additionally for each method, the number of iterations used is provided in Table 1. Table 1

AN US

Some details of the five SR algorithms investigated and our proposed method SR method

Iteration number

predetermined parameters

Tuning

VTV[21]

Step1 Step2 20 15 10 50

α=0.9

regularized parameters ---regularized parameters regularized parameters

RSATV[13]

M

α=0.8,q=2, c=1, α=0.8,β=0.05,τ=20 Noise variance=σ2

----

̂

bi-cubic interpolation of the first frame

̂ For all

CE

Proposed method

80(srf=2) 120(srf=4) 80(srf=2) 120(srf=4) Automatic determination

ED

AN+MBTV[12]

σ2=10,ε2=0.5 q=2

PT

KFB[23] L1+BTV[11]

Initial values

AC

4-1- Results of the simulated frames 4-1-1 Effectiveness of the proposed method in the fine details restoring In this section, to validate the ability of the proposed method in details restoring, an experiment is designed conducted to a test pattern as the HR image containing different frequencies in spatial domain (Fig.3(a)). The process of the observation model is carried out in the following steps. At first, the original HR image is shifted in horizontal and vertical 12

ACCEPTED MANUSCRIPT planes while the translations are chosen randomly from the ranges between (0-1). Then, the shifted images are blurred with a low pass filter as imaging system PSF listed in Table 2 (srf=2). In the third step, the blurred image is down sampled by factor 2. Finally, zero mean Gaussian noise with variance σ2=4 is added to obtain the LR image. In this experiment, 12

Table 2 Characteristics of three types of PSFs used in our experiments

srf=4

size Strength(∑

)

size Strength(∑

PSF2 Gaussian 3×3, σ=1 0.1256 5×5, σ=1 0.8250

AN US

srf=2

PSF1 Square 2×2 0.2500 4×4 0.6250

CR IP T

LR images are provided.

)

PSF3 Gaussian 3×3, σ=2 0.1119 5×5, σ=2 0.4330

The experimental quantity results are reported in Table 3 including PSNR, SSIM [25] and

M

FSIM [26] of VTV [21], KFB [23], L1+BTV [11], AN+MBTV [12], and RSATV [13] methods. Added to which is our proposed method with different PSFs. Moreover, some

ED

reconstructed HR images are given in Figs. 3 through 5.

PT

It can be seen that for the cropped regions, the proposed method achieves higher PSNR, SSIM and FSIM than that of other methods. In the other hand, the results obtained by using

CE

the proposed method are more appealing in preserving the fine details containing thin lines and corners in the test pattern. Moreover, for the total reconstructed test pattern, our proposed

AC

method restores higher quality images based on SSIM and FSIM. Fig.4 and 5 show the corresponding visual effects.

13

ACCEPTED MANUSCRIPT Table 3 PSNR, SSIM, and FSIM of reconstructed test pattern (cropped region) by VTV [21], KFB [23], L1+BTV [11], AN+MBTV[12], RSATV[13], and our proposed methods with different PSFs Assessment

VTV [21]

L1+BTV [11]

KFB [23]

AN+MBTV[12]

Proposed

RSATV[13]

20.54(12.70)

26.53(17.78)

19.32(15.89)

30.23(21.59)

29.85(23.98)

30.80(29.62)

SSIM

0.9792(0.7096)

0.9732(0.9042)

0.9066(0.8793)

0.9748(0.9663)

0.9934(0.9875)

0.9978(0.9902)

FSIM

0.9182(0.7418)

0.9869(0.9541)

0.9636(0.8527)

0.9910(0.9813)

0.9973(0.9757)

0.9999(0.9969)

PSNR

17.08(7.93)

24.18(13.66)

--

22.64(12.70)

24.05(12.39)

22.41(16.81)

SSIM

0.9313(0.0785)

0.9676(0.7546)

---

0.9651(0.6871)

0.9651(0.6833)

0.9790(09117)

FSIM

0.8956-0.2446

0.9824-0.8016

---

0.9795-0.7365

0.9844-0.7547

0.9945-0.9035

PSNR

16.42(7.83)

21.57(10.15)

---

21.91(10.86)

22.84(10.83)

21.64(14.04)

SSIM

0.9127(0.0605)

0.9481(0.4697)

---

0.9732(0.5476)

0.9535(0.5729)

0.9538(0.7117)

FSIM

0.8945(0.2325)

0.9727(0.4470)

---

0.9811(0.5081)

0.9769(0.4936)

0.9853(0.8182)

ED

M

AN US

CR IP T

PSNR

(b)

(c)

(e)

(f)

CE

PT

(a)

AC

PSF3

PSF2

PSF1

method

(d)

Fig.3, Experimental results of the test pattern image by using PSF1: (a) Original HR image. Reconstructed images of the: (b) VTV[21], (c) L1+BTV[11],(d) AN+MBTV[12], (e) RSATV[13], and (f) The proposed SR methods.

14

ACCEPTED MANUSCRIPT

(b)

(d)

(e)

(c)

CR IP T

(a)

(f)

(b)

(c)

(d)

(e)

(c)

(d)

(e)

CE

(a)

PT

ED

(a)

M

AN US

Fig.4, Detailed region cropped from Fig.3.

(b)

AC

Fig.5, Detailed region cropped from reconstructed HR images obtained by: (a) VTV [21] ,(b) KFB[23],(c) L1+ BTV[11],(d) AN+MBTV[12], (e) RSATV[13], and (f) The proposed SR methods. First row and second row corresponding to PSF2 and PSF3, respectively.

4-1-2- Testing our proposed method conducted to traditional images In this part of simulation experiments, three popular images are chosen, Boat, Zebra, and Barbara of a 400 400 size as HR images having proper details. The process of the observation model is carried out in the same manner which is described in the previous

15

ACCEPTED MANUSCRIPT section in the cases where srf=2. For the cases where srf=4, the translations are chosen randomly from ranges between (0-4), imaging system PSF is selected from the second row of the Table 2, and the blurred image is down sampled by factor 4. 30 LR images are generated in the cases where srf=4.

CR IP T

The experimental quantity results are reported in Table 4 and 5 including PSNR, SSIM [25] and FSIM [26] of VTV [21], KFB [23], L1+BTV [11], AN+MBTV [12],and RSATV [13] methods. Added to which is our proposed method with different PSFs where the SR factors (srf) are considered 2 and 4. Moreover, some reconstructed HR images, while cropped

AN US

regions, and corresponding metric indexes are given in Figs. 6 through 8 in the case where srf=2 and, in Figs. 9 through 11 in the case where srf=4, respectively. The following points deserve attention:

Based on the three metrics PSNR, SSIM, and FSIM, for most cases, the obtained

M



images by our proposed method exhibit higher quality measures computed with those

The edges, the characters, and the lines in the reconstructed images accomplished via

PT



ED

of other methods, regardless of the srf, the image, and the PSF.

the proposed method are much more vivid and distinct compared with those of other

CE

methods. These are reflected in the shown cropped regions and corresponding metrics.

The output quality of the VTV [21] method is completely damaged in case where

AC



srf=4.



For the image containing more details such as Zebra and Barbara, the proposed method is more effective i.e. the details are better reconstructed through our method.

16

ACCEPTED MANUSCRIPT 

Based on the obtained results of the quantity and visual dimensions, our reconstructed output quality- corresponding to PSF1and PSF2 in the cases where srf=2 and srf=4 respectively- surpasses all other considered cases.



The order of the outputs qualities corresponding to the PSFs are as follows: PSF1>PSF2>PSF3 and PSF2>PSF1>PSF3 in the cases where srf=2 and srf=4,

CR IP T

respectively which are compatible with strength ordering of the PSF filters. It is noted that the strength of the PSF filters (∑

) are presented in Table 2.

In general, based on the reported metrics in Table 4 ,Table 5 and the images displayed in

AN US

Figs. 6 through 11, comparatively, the reconstructed HR images and the restored details produced by our proposed method are superior to those obtained through the five other methods in cases where srf=2 and srf=4.

M

In addition, one of the most important advantages of the proposed method is that it does not require any parameter tuning, while generally speaking, in the case of the regularized SR

ED

methods like the ones utilized in our experiments, it is necessary to adjust the tuning parameters-dependent upon the image properties-effected by several rounds of running until

PT

the acceptable output quality is attained. Moreover, the KFB [23] method is limited by the

AC

CE

size of PSF while our proposed method supports PSFs of whatever size and form.

17

ACCEPTED MANUSCRIPT Table 4 PSNR, SSIM, and FSIM of reconstructed HR images by VTV [21], KFB [23], L1+BTV [11], AN+MBTV[12], RSATV[13], and our proposed methods with different PSFs in the case where srf=2 Zebra

Barbara

PSF3

PSF1

PSF2

PSF3

PSF1

PSF2

PSF3

VTV [21]

34.23

32.38

31.51

30.62

27.94

26.51

31.23

27.14

26.00

KFB[ 23]

33.82

---

---

30.46

---

---

31.50

---

---

L1+BTV [11]

34.07

32.34

31.99

30.53

28.09

27.08

31.04

27.73

26.94

AN+MBTV [12]

36.32

33.87

33.88

33.47

29.64

RSATV [13]

35.44

33.02

32.87

32.47

29.21

Proposed Method

35.92

33.11

33.00

35.27

29.99

CR IP T

PSF2

28.70

33.69

29.08

28.58

28.65

32.23

28.76

28.29

29.04

35.74

29.15

28.65

VTV [21]

0.9952 0.9890 0.9835 0.9956 0.9858 0.9753 0.9898 0.9542 0.9268

KFB [23]

0.9930

---

---

0.9888

---

---

0.9882

---

---

0.9873 0.9847 0.9806 0.9877 0.9767 0.9708 0.9781 0.9523 0.9321

AN+MBTV [12]

0.9957 0.9923 0.9902 0.9969 0.9920 0.9851 0.9927 0.9757 0.9621

RSATV [13]

0.9937 0.9877 0.9844 0.9958 0.9907 0.9856 0.9906 0.9698 0.9555

Proposed Method

0.9960 0.9905 0.9866 0.9986 0.9940 0.9879 0.9958 0.9781 0.9616

VTV[21]

0.9967 0.9927 0.9914 0.9959 0.9892 0.9846 0.9943 0.9787 0.9723

KFB [23]

0.9950

---

AN US

L1+BTV [11]

---

0.9908

---

---

0.9928

---

---

0.9913 0.9907 0.9893 0.9914 0.9839 0.9822 0.9880 0.9758 0.9693

AN+MBTV [12]

0.9980 0.9962 0.9954 0.9970 0.9929 0.9916 0.9960 0.9884 0.9841

RSATV [13]]

0.9972 0.9946 0.9932 0.9951 0.9904 0.9895 0.9953 0.9864 0.9817

Proposed Method

0.9986 0.9959 0.9940 0.9989 0.9941 0.9878 0.9980 0.9885 0.9806

ED

M

L1+BTV [11]

AC

CE

PT

FSIM

SSIM

PSNR

Boat PSF1

18

ACCEPTED MANUSCRIPT

PSNR=32.52 SSIM=0.9486 FSIM=0.9764

PSNR=32.34 SSIM=0.9430 FSIM=0.9776

(b)

M

(c)

PSNR=35.03, SSIM=0.9663 FSIM=0.9869

AN US

PSNR=32.52 SSIM=0.9364 FSIM=0.9772

CR IP T

(a)

(e)

ED

PSNR=34.22, SSIM=0.9606 FSIM=0.9846

(d)

PSNR=35.32, SSIM=0.9516 FSIM=0.9914

(f)

AC

CE

PT

Fig.6, Reconstructed images , regional details, and corresponding metric indexes of Zebra by using: (a) VTV [21] ,(b) KFB[23], (c) L1+BTV[11],(d) AN+MBTV[12], (e) RSATV[13], and (f) The proposed SR methods while considering PSF1 in the case where srf=2.

19

ACCEPTED MANUSCRIPT

PSNR=23.86, SSIM=0.4420, FSIM=0.7027

PSNR=24.53, SSIM=0.5169, FSIM=0.7490

(b)

(d)

ED

M

(c)

PSNR=25.29, SSIM=0.5643, FSIM=0.7993

AN US

PSNR=25.19, SSIM=0.5526, FSIM=0.7759

CR IP T

(a)

PSNR=25.71, SSIM=0.6040, FSIM=0.8031

(e)

AC

CE

PT

Fig.7, Reconstructed images , regional details, and corresponding metric indexes of Barbara by using: (a) VTV [21] , (b) L1+BTV[11],(c) AN+MBTV[12], (d) RSATV[13], and (e) The proposed SR methods while considering PSF2 in the case where srf=2.

20

ACCEPTED MANUSCRIPT

PSNR=28.31, SSIM=0.7903, FSIM=0.8778

PSNR=23.86, SSIM=0.4420, FSIM=0.7027

(b)

CR IP T

(a)

PSNR=29.24, SSIM=0.8308, FSIM=0.9053

AN US

PSNR=29.93, SSIM=0.8434, FSIM=0.9009

(d)

ED

M

(c)

PSNR=30.44, SSIM=0.8565, FSIM=0.9158

(e)

AC

CE

PT

Fig.8, Reconstructed images, regional details, and corresponding metric indexes of Boat by using: (a) VTV [21], (b) L1+BTV[11], (c) AN+MBTV[12], (d) RSATV[13], and (e) The proposed SR methods while considering PSF3 in the case where srf=2.

21

ACCEPTED MANUSCRIPT Table 5 PSNR, SSIM, and FSIM of reconstructed HR images by VTV [21], KFB [23], L1+BTV [11], AN+MBTV [12], RSATV [13], and our proposed methods with different PSFs in the case where srf=4 Boat Zebra Barbara PSF3

PSF1

PSF2

PSF3

PSF1

PSF2

PSF3

29.185

---

---

23.807

---

---

24.973

---

---

29.04

29.48

28.21

23.58

24.73

22.72

24.61

25.25

23.81

AN+MBTV [12] 31.67

32.45

30.26

26.30

28.19

24.78

26.15

27.19

24.99

31.14

31.79

29.92

25.77

27.51

24.24

25.99

26.33

24.81

Proposed Method 30.30

31.013

29.21

25.58

27.76

24.08

25.79

27.35

24.62

CR IP T

KFB [ 23] L1+BTV [11] RSATV [13]

SSIM

PSF2

25.182 25.340 25.476 20.452 20.860 20.371 22.910 22.810 22.808

VTV [21]

0.8685 0.8859 0.8919 0.7906 0.7744 0.7897 0.7980 0.7915 0.7992

KFB [23]

0.9558

L1+BTV [11]

---

---

0.9375

---

---

0.8929

---

---

0.9542 0.9618 0.9450 0.9188 0.9521 0.8884 0.8780 0.9084 0.8478

AN+MBTV [12] 0.9768 0.9864 0.9646 0.9591 0.9840 0.9373 0.9181 0.9477 0.8828 RSATV [13]

0.9678 0.9798 0.9558 0.9541 0.9811 0.9307 0.9048 0.9353 0.8717

AN US

PSNR

PSF1 VTV [21]

FSIM

Proposed Method 0.9639 0.9591 0.9488 0.9529 0.9838 0.9340 0.9072 0.9465 0.8741 VTV [21]

0.9222 0.9303 0.9378 0..8863 0.9012 0.8956 0.9309 0.9176 0.9296

KFB [23]

0.9843

L1+BTV [11]

---

---

0.9788

---

---

0.9732

---

---

0.9807 0.9798 0.9763 0.9720 0.9713 0.9579 0.9653 0.9596 0.9555

AN+MBTV [12] 0.9933 0.9940 0.9882 0.9902 0.9918 0.9812 0.9812 0.9792 0.9688 RSATV [13]

0.9917 0.9925 0.9861 0.9865 0.9887 0.9750 0.9775 0.9762 0.9656

AC

CE

PT

ED

M

Proposed Method 0.9905 0.9924 0.9856 0.9876 0.9902 0.9777 0.9790 0.9792 0.9668

22

ACCEPTED MANUSCRIPT

PSNR=20.12, SSIM=0.2058, FSIM=0.5666

PSNR=20.42, SSIM=0.2543, FSIM=0.6112

(b)

(d)

ED

M

(c)

PSNR=20.75, SSIM=0.3100, FSIM=0.6456

AN US

PSNR=20.13, SSIM=0.2201, FSIM=0.5628

CR IP T

(a)

PSNR=20.81, SSIM=0.3274, FSIM=0.6806

(e)

PSNR=21.16, SSIM=0.3689, FSIM=0.6872

(f)

AC

CE

PT

Fig.9, Reconstructed images , regional details, and corresponding metric indexes of Barbara by using: (a) VTV [21] ,(b) KFB[23], (c) L1+BTV[11],(d) AN+MBTV[12], (e) RSATV[13], and (f) The proposed SR methods while considering PSF1 in the case where srf=4.

23

ACCEPTED MANUSCRIPT

PSNR=24.00, SSIM=0.6855, FSIM=0.8080

PSNR=26.74, SSIM=0.7706, FSIM=0.8487

(b)

CR IP T

(a)

PSNR=28.30, SSIM=0.7955, FSIM=0.8742

AN US

PSNR=28.54, SSIM=0.8087, FSIM=0.8740

(d)

ED

M

(c)

PSNR=28.55, SSIM=0.7771, FSIM=0.8905

(e)

AC

CE

PT

Fig.10, Reconstructed images , regional details, and corresponding metric indexes of Zebra by using: (a) VTV [21] , (b) L1+BTV[11],(c) AN+MBTV[12], (d) RSATV[13], and (e) The proposed SR methods while considering PSF2 in the case where srf=4.

24

ACCEPTED MANUSCRIPT

PSNR=25.26, SSIM=0.6837, FSIM=0.8176

PSNR=22.79, SSIM=0.5848, FSIM=0.7794

(b)

CR IP T

(a)

PSNR=26.98, SSIM=0.7617, FSIM=0.8718

AN US

PSNR=26.67, SSIM=0.7510, FSIM=0.8581

(d)

ED

M

(c)

PSNR=30.44, SSIM=0.8565, FSIM=0.9158

(e)

CE

PT

Fig.11, Reconstructed images , regional details, and corresponding metric indexes of Boat by using: (a) VTV [21] , (b) L1+BTV[11], (c) AN+MBTV[12], (d) RSATV[13], and (e) The proposed SR methods while considering PSF3 in the case where srf=4.

AC

4-2- Testing the proposed method under different noise intensities In order to evaluate the effectiveness of our proposed SR method under different conditions of noise intensity, we carry out experiments like the ones we did in the previous section except for the fact that this time different noise values are added. Firstly, simulated datasets of Zebra and Barbara images are generated while adding noise under different variances 4, 10, 16, 20 and 28. Then we apply the SR methods including our proposed method and two SR methods AN+MBTV [12] and RSATV[13] in the case where srf=2 using PSF1. The 25

ACCEPTED MANUSCRIPT results obtained from the quantitative evaluation results for reconstructed HR images and cropped regions are shown in Tables 6 and 7 while some images are shown in Figs.12 and 13. It can be deduced that: 

In most cases, based on the three reported metrics, the proposed method generates a



CR IP T

higher quality HR image than that of the other ones. When there is an increase in noise, in all cases, the three index of the cropped regions of the HR images -obtained through the proposed method- is higher than that of other methods even in the case of AN+MBTV [12] method. This means that the features

AN US

and the details of the HR images are better restored through our proposed method, a fact reflected by the FSIM metric in the high noisy cases.

PSNR SSIM FSIM PSNR SSIM FSIM PSNR SSIM FSIM PSNR SSIM FSIM PSNR SSIM FSIM

PT

4

Index

10

CE

16

AC

20

28

AN+MBTV[12]

ED

Noise Variance

M

Table 6 Quantitative evaluation results of different noise intensities for Zebra image (cropped region) experiment

33.66(32.97) 0.9972(0.9087) 0.9972(0.9476) 33.04(31.99) 0.9959(0.8881) 0.9964(0.9381) 32.57(31.57) 0.9946(0.8823) 0.9957(0.9317) 32.29(30.96) 0.9936(0.8679) 0.9951(0.9256) 31.79(31.11) 0.9921(0.8677) 0.9943(0.9284)

26

RSATV[13]

Proposed method

32.47(32.30) 0.9958(0.8964) 0.9951(0.9478) 31.84(31.68) 0.9938(0.8848) 0.9937(0.9440) 31.35(30.86) 0.9923(0.8683) 0.9925(0.9347) 31.11(30.52) 0.9913(0.8552) 0.9921(0.9326) 30.65(30.25) 0.9891(0.8536) 0.9906(0.9285)

35.19(34.30) 0.9982(0.9206) 0.9986(0.9569) 33.19(32.87) 0.9959(0.8969) 0.9969(0.9502) 32.08(31.75) 0.9937(0.8737) 0.9955(0.9380) 31.54(30.99) 0.9926(0.8491) 0.9946(0.9324) 30.77(30.40) 0.9899(0.8479) 0.9931(0.9289)

ACCEPTED MANUSCRIPT Table 7 Quantitative evaluation results of different noise intensities for Barbara image (cropped region) experiment

10

16

20

RSATV[13]

Proposed method

PSNR SSIM FSIM PSNR SSIM FSIM PSNR SSIM FSIM PSNR SSIM FSIM PSNR SSIM FSIM

33.69(31.20) 0.9927(0.9382) 0.9960(0.9605 33.10(30.47) 0.9898(0.9340) 0.9948(0.9547) 32.50(30.38) 09857(0.9288 0.9929(0.9539) 32.03(29.66) 0.9832(0.9065) 0.9915(0.9416) 31.56(29.38) 0.9793(0.9018) 0.9889(0.9389)

33.24(30.59) 0.9906(0.9412) 0.9953(0.9614) 32.35(30.47) 0.9839(0.9363) 0.9924(0.9581) 31.74(30.36) 0.9775(0.9320) 0.9899(0.9564) 31.32(30.06) 0.9735(0.9271) 0.9882(0.9539) 30.85(29.87) 0.9669(0.9175) 0.9855(0.9475)

35.70(35.31) 0.9947(0.9725) 0.9975(0.9776) 33.50(33.08) 0.9881(0.9572) 0.9946(0.9684) 32.33(32.08) 0.9825(0.9481) 0.9920(0.9619) 31.83(31.16) 0.9791(0.9382) 0.9906(0.9555) 31.00(30.42) 0.9724(0.9247) 0.9877(0.9500)

AC

CE

PT

ED

M

28

AN+MBTV[12]

CR IP T

4

Index

AN US

Noise Variance

(a)

(b)

(c)

Fig.12, One of the generated LR sequences under different noise intensities. Noise variance equals to: (a) 4, (b) 16, and (c) 28. 27

AN US

CR IP T

ACCEPTED MANUSCRIPT

(a)

(b)

(c)

Fig.13, Experimental results of the Zebra image by using PSF1 in the case where srf=2 and noise variance=16. Reconstructed images (first row) and regional details (second row) of: (a) AN+MBTV[12], (b) RSATV[13], and (c) The proposed SR methods.

M

4-3- Effects of the initial and the setting values on the performance of the proposed method

ED

The Kalman filter developed through Eqs.(11), (12), and (13) requires only one setting value

considered as

PT

for the noise variance and one initial value for the covariance matrix, which is generally where

and

are the identity matrix and the

CE

initial scalar value respectively. First, to analyze the sensitivity of the performance to parameter

, in both cases where srf=2 using PSF1 and where srf=4 using PSF2, for Barbara

AC

datasets, the algorithm is repeated with different values of; subsequently, the obtained PSNRs versus the iteration schemes are plotted in Fig.14(a). It can be observed that by considering a lower value for

, only the required number of iterations increases while the algorithm

converges to the same final PSNR in all cases. This confirms that the robustness of our proposed algorithm to initialization. Based on the obtained results, in all experiments and

are considered in the cases where srf=2 and srf=4, respectively.

28

ACCEPTED MANUSCRIPT Second, to verify the robustness of the performance of proposed SR method to setting noise variance, the algorithm is repeated as in the previous paragraph except that the setting noise variance is considered differently. Subsequently the obtained PSNRs versus iteration number are plotted in Fig.14(b). It can be seen that when setting noise variance is changed, only the

AN US

CR IP T

required iteration number is changed while the final PSNR remains constant.

CE

PT

ED

M

(a)

(b)

Fig.14, PSNR of reconstructed Barbara image versus iteration number considering different:

AC

(a) Initial values for

: (b) Setting values for noise variance. Left and right columns

corresponding to cases where srf=2 using PSF1 and srf=4 using PSF2, respectively. 4-4-Real-world frames To verify the effectiveness of our algorithm in practical applications, we also conducted some experiments with two real-world images, including Alpaca and Emily sequences [27]. These sequences have adequate details which can visually differentiate the results of our proposed

29

ACCEPTED MANUSCRIPT method against the five other methods compared. At the beginning, 30 frames of each sequence (16-45 from Alpaca and 26-55 from Emily) are selected as the LR with a specified size of 96 128. The well-known block-based registration algorithm is also used to estimate the relevant translations with sub-pixel accuracy [28].

CR IP T

Some reconstructed HR images and their relevant zoomed regional details are also provided in Figs. 15 through 22 in the case where srf=2 and in Figs. 23 through 30 in the case where srf=4, respectively. To exactly show the restored details, we crop two regions for every reconstructed HR image containing middle and finer details. With regard to Figs. 15 through



AN US

30, the following can be deduced:

In the case where srf=2, the details are clearly restored by our proposed method as far as PSF1 is concerned. To take an example, two characters 'i' and 'e' are noticeable in



M

Fig. 16.

The characters and the numbers in the zoomed regional details restored by

ED

AN+MBTV [12],RSATV [13] and our proposed method are very close to each other in the case where srf=2 when PSF2 and PSF3 are selected. As in the case where srf=4, when PSF1 is used, the reconstructed HR image restored

PT



CE

by RSATV [13] and our proposed method are very close to each other and better than that produced via other methods. The restored lines, the edges, characters, numbers and other middle details by

AC



AN+MBTV[12] ,RSATV[13] and our proposed method are clear with minimum dissimilarity as the two connect for the case where srf=4 and PSF2 is considered.

Altogether, the middle and the finer details of the HR images are better restored in respect of resolution by our proposed method when PSF1 and PSF2 are selected in cases where srf=2 and srf=4, respectively. (cf. Figs. 16, 22, 26, and 30). Further, these cases are compatible with 30

ACCEPTED MANUSCRIPT the simulation results acquired in the previous section. This means that the SR reconstruction

(b)

(c)

AN US

(a)

CR IP T

performance of our proposed method is far superior.

(d)

(e)

(f)

AC

CE

PT

ED

M

Fig.15, Experimental results of the real-world Alpaca images using PSF1. Reconstructed HR image by: (a) VTV [21], (b) KFB[23], (c) L1+BTV[11],(d) AN+MBTV[12], (e) RSATV[13], and (f) The proposed methods in the case where srf=2.

(a)

(b)

(c)

(d)

(e)

(f)

Fig.16, Two sets of detailed regions cropped from Fig.15. First row Alp-Region1, second row Alp-Region2.

31

ACCEPTED MANUSCRIPT

(a)

(c)

AN US

CR IP T

(b)

(d)

(e)

AC

CE

PT

ED

M

Fig.17, Experimental results of the real-world Alpaca images using PSF2. Reconstructed HR image by: (a) VTV[21], (b) L1+BTV[11], (c) AN+MBTV[12], (d) RSATV[13], and (e) The proposed methods in the case where srf=2.

(a)

(b)

(c)

(d)

(e)

Fig.18, Two sets of detailed regions cropped from Fig.17. First row Alp-Region1, second row Alp-Region2.

32

ACCEPTED MANUSCRIPT

(b)

(c)

CR IP T

(a)

(e)

(d)

AC

CE

PT

ED

M

AN US

Fig.19, Experimental results of the real-world Alpaca images using PSF3. Reconstructed HR image by: (a) VTV[21],(b) L1+BTV[11],(c) AN+MBTV[12], (d) RSATV[13], and (e) The proposed methods in the case where srf=2.

(a)

(b)

(c)

(d)

(e)

Fig.20, Two sets of detailed regions cropped from Fig.19. First row Alp-Region1, second row Alp-Region2.

33

ACCEPTED MANUSCRIPT

(b)

(d)

(e)

(c)

CR IP T

(a)

(f)

CE

PT

ED

M

AN US

Fig.21, Experimental results of the real-world Emily images using PSF1. Reconstructed HR image by: (a) VTV [21], (b) KFB[23], (c) L1+BTV[11](d) AN+MBTV[12], (e) RSATV[13], and (f) The proposed methods in the case where srf=2.

(a)

(b)

(d)

(c)

(e)

(f)

AC

Fig.22, Two sets of detailed regions cropped from Fig.21. First row Emil-Region1, second row Emil-Region2.

34

ACCEPTED MANUSCRIPT

(b)

(d)

(e)

(c)

CR IP T

(a)

(f)

AC

CE

PT

ED

M

AN US

Fig.23, Experimental results of the real-world Alpaca images using PSF1.Reconstructed HR image by: (a) VTV[21] , (b) KFB[23], (c) L1+BTV[11], (d) AN+MBTV[12], (e) RSATV[13], and(f) The proposed methods in the case where srf=4.

(a)

(b)

(d)

(c)

(e)

(f)

Fig.24, Two sets of detailed regions cropped from Fig.23. First row Alp-Region1, second row Alp-Region2.

35

ACCEPTED MANUSCRIPT

(a)

(c)

CR IP T

(b)

(d)

(e)

AC

CE

PT

ED

M

AN US

Fig.25, Experimental results of the real-world Alpaca images using PSF2. Reconstructed HR image by: (a) VTV[21], (b) L1+BTV[11],(c) AN+MBTV[12], (d) RSATV[13] , and (e) The proposed methods in the case where srf=4.

(a)

(b)

(c)

(d)

(e)

Fig.26, Two sets of detailed regions cropped from Fig.25. First row Alp-Region1, second row Alp-Region2.

36

ACCEPTED MANUSCRIPT

(a)

(c)

CR IP T

(b)

(d)

(e)

CE

PT

ED

M

AN US

Fig.27, Experimental results of the real-world Alpaca images using PSF3. Reconstructed HR image by: (a) VTV [21], (b) L1+BTV[11], (c) AN+MBTV[12], (d) RSATV[13], and (e) The proposed methods in the case where srf=4.

(a)

(b)

(c)

(d)

(e)

AC

Fig.28, Two sets of detailed regions cropped from Fig.27. First row Alp-Region1, second row Alp-Region2.

37

ACCEPTED MANUSCRIPT

(b)

(c)

CR IP T

(a)

(d)

(e)

AC

CE

PT

ED

M

AN US

Fig.29, Experimental results of the real-world Emily images using PSF2.Reconstructed HR image by: (a) VTV[21], (b) L1+BTV[11], (c) AN+MBTV[12], (d) RSATV[13], and (e) The proposed methods in the case where srf=4.

(a)

(b)

(c)

(d)

(e)

Fig.30, Two sets of detailed regions cropped from Fig.29. First row Emil-Region1, second row Emil-Region2.

4-5-Comparing the running time of different methods In order to fairly compare the computational complexity of our proposed method with that of other five SR methods, the measured running times for a real sequence Alpaca in two cases

38

ACCEPTED MANUSCRIPT where srf=2 and srf=4 are summarized in Table 8. It should be noted that the experiments are implemented under Matlab R2016 environment, exploiting Core-i7 2.2 GHz CPU. It is seen that: 

Our proposed method converges toward the best results in the reasonable running



CR IP T

time being lower than that of the KFB[23] method. Although the proposed method displays higher running time in comparison with the VTV[21], L1+BTV[11], and RSATV[13] methods, based on the previous results obtained and reported in Tables 3 through 7,our proposed method achieves higher

AN US

quality HR images.

It is noted that our proposed method can be accelerated by parallel implementation process

M

while the other methods are unfeasible.

ED

Table 8 Execution time (in seconds) of our proposed and other five SR methods

4.25

94.123

15.87

52.12

27.29

54.31

CE

srf=4

PT

Case VTV[21] KFB[23] L1+BTV[11] AN+MBTV[12] RSATV[13] Proposed method srf=2 0.63 50.40 3.746 13.74 6.85 28.11

AC

4-6- Scalable SR by partial spectrum rejection technique In some of the SR applications, the desired HR image size is quite fixed, that is, impossible to change. In other words, the size of the required HR image is predetermined by the resolution and the aspect ratio of the monitoring system. Employing the proposed SR algorithm, we can reconstruct the HR image of any size

,

in a very simple manner, while

other traditional SR methods are unfeasible with regard to the HR image size and scale. In

39

ACCEPTED MANUSCRIPT this respect, considering our proposed algorithm, it is sufficient to reject some of the end spectrum of ̂ ̂

as follows (cf. Appendix B): ̂

(15)

;

Fig.31 illuminates how the end of spectrum of the HR image implements the rejection

M

AN US

CR IP T

process before using inverse DFT.

ED

Fig.31, End of spectrum rejection to obtain HR image with an arbitrary size and scale.

In general, in order to obtain HR images with any size and scale, a well-known technique is

PT

to conduct one of the traditional resampling methods after SR procedure. In order to evaluate

CE

the proposed scalable SR method, two traditional resampling methods employed executed by imresize function in Matlab and an HR image the same size as that of the proposed method is

AC

obtained.Figs.32 and 33 present the reconstructed HR images of the two previous real world sequences of three different sizes and scales obtained via both partial spectrum rejection technique and resampling method which drawn upon nearest and bicubic interpolation. As is shown in the two figures, the images reconstructed by the scalable SR method exhibit lower noises because of rejecting high frequency components containing major noise elements. In addition, the time is saved through our proposed technique since no resampling process is resorted. 40

ACCEPTED MANUSCRIPT Size: 370 500 (real scale)

Size: 250 500

PT

ED

M

AN US

CR IP T

Size: 350 350

AC

CE

Fig. 32, Reconstructed images of Alpaca displaying different sizes and scales using partial spectrum rejection technique (first row), resampling within nearest interpolation method (second row) ,and resampling within bicubic interpolation(third row) .

41

ACCEPTED MANUSCRIPT Size: 370 500 (real scale)

Size: 250 500

ED

M

AN US

CR IP T

Size: 350 350

CE

PT

Fig.33, Reconstructed images of Emily displaying different sizes and scales using partial spectrum rejection technique (first row), resampling within nearest interpolation method (second row) , and resampling within bicubic interpolation(third row) .

AC

5- Conclusion

Applying the Kalman filter to the frequency domain, we have introduced a multi-frame super-resolution algorithm. Unlike previous studies, exploiting similar methods, we have considered a method which solves two major problems of the conventional Kalman filtering technique requiring neither approximations nor limitations on space-invariant PSF of any type, whatsoever. Further, parallel implementation enables us to speed up our proposed

42

ACCEPTED MANUSCRIPT method making it feasible to set PSF of any form and size, which is a great help in attaining higher quality HR images. The presented technique requires no numerical manual parameters or tuning. Experiments are conducted on both simulated and real data. Because the proposed method emphasizes restoring details which is the main task of all SR techniques, it is more effective to reconstruct images containing finer particulars. The obtained results are also a

CR IP T

confirmation of the fact that the algorithm offered here outperforms those developed by other approaches. Additionally, this algorithm promises better results as regards quantity and vision compared with those of similar methods. Further, a scalable SR algorithm is advanced through partial spectrum rejection technique, as a result of which an HR image of an arbitrary

AN US

size and scale is obtained.

Although, the SR method set forth in this research is applicable to a wide swath of real data, developing the latter method so that it incorporates such more complicated motion conditions

ED

M

as rotation and zoom is a part of our future research project.

PT

Appendix A

To simplify, we consider one dimensional signal

CE

Gaussian noise and ∑

, as a white

as its FT in the following relation:

(A1)

It is not difficult to show that

is a zero mean Gaussian noise. Now, let us calculate the

AC

,

relevant autocorrelation:

̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ∑

∑

̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ∑ ∑

̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅

(

)

(A2) 43

ACCEPTED MANUSCRIPT

∑

(

∑

Where * is a conjugate,

)

(

∑

represents power of noise

function. Relation A2 shows that

)

and

stand for impulse

is a white Gaussian noise.

Appendix B

is assumed to be sampled at or above Nyquist rate.

)

,

(B1)

( Where

)

AN US

(

( ) which

CR IP T

( ) is a one dimensional continuous band limited signal with FT

Supposing

,

is the sampling period. Now, let us cut a part of spectrum of

as follows:

(B3)

M

,

(B2)

its corresponding signal is obtained by inverse DFT as is borne out by the

ED

Where

∑

(

)

(B4)

AC

CE

∑

PT

following:

References

[1] L. Yue, H. Shen, L. Li, Q. Yuan, H. Zhang, Image super-resolution: the techniques, application, and future, Signal Process. 128 (2016) 389-408.

[2] A. Rahimi, P. Moallem, K. Shahtalebi, M. Momeni, Preserving quality in minimum frame selection within multi-frame super-resolution, Digital Signal Process. 72 (2018) 19-43. [3] T. S. Huang, R. Y. Tsai, Multi-frame image restoration and registration, Advances in Computer Vision and Image Process. 1, JAI Press Inc., Greenwich, CT (1984) 317–339.

44

ACCEPTED MANUSCRIPT [4] A. M. Tekalp, M. K. Ozka, and M. I. Sezan, High-resolution image reconstruction from lowerresolution image sequences and space-varying image restoration, in Proc. IEEE Int. Conf. Acoustics, Speech, and Signal Processing, San Francisco, California, 3(1992)169-172. [5] S. P. Kim, N. K. Bose, and H. M Valenzuela, Recursive reconstruction of high resolution image from noisy under sampled multiframes, IEEE Trans. on Speech and Signal Process. 38 (1990) 1013–027.

reconstruction algorithm, Opt. Eng. 24 (2003) 1408-1432.

CR IP T

[6] S. Rhee, M G. Kang, Discrete cosine transform based regularized high-resolution image

[7] N. Nguyen, P. Milanfar, G. Colub, A computationally efficient superresolution image reconstruction algorithm, IEEE Trans. Image Process. 10(2001) 573-583.

[8] M. Irani, S. Peleg, Improving resolution by image registration, CVGIP: Graph. Model, Image

AN US

Process. 53(1991) 231-239.

[9] A.M. Tekalp, M.K. Ozkan, M.I. Sezan, High-resolution image reconstruction from lowerresolution image sequences and space-varying image restoration, in: Proceeding of the IEEE International Conference on Acoustics, Speech, and Signal Processing. 1992. pp. 169-172. [10] H. Shen, I. ZHANG, B. Huang, et al., A MAP approach for joint motion estimation,

M

segmentation, and super-resolution IEEE Trans. Image Process. 16(2007) 479-490. [11] S. Farsiu, M.D. Robinson, M. Elad, P. Milanfar, Fast and robust multi-frame super resolution,

ED

IEEE Trans. Image Process. 13 (2004) 1327–1344. [12] X. Zeng, L. Yang, A robust multiframe super-resolution algorithm based on half-quadratic

PT

estimation with modified BTV regularization, Digital Signal Process. 23 (2013) 98-109. [13] Q. Yuan, L. Zhang, H. Shen, and P. Li, Regional spatially adaptive total variation super-

CE

resolution with spatial information filtering and clustering, IEEE Trans. Image process. 22 (2013) 2327-2342.

AC

[14] L. Yue, H. Shen, Q. Yuan, and L. Zhang, A locally adaptive L1-L2 norm for multi-frame superresolution of images with mixed noise and outliers, Signal process. 105 (2014) 156-174.

[15] X. Yang, T. Liu,D. Zhou, A multi-frame adaptive super-resolution method using double channel and regional pixel information, Optik 126 (2015) 5850-5858.

[16] H. Shen, L. Peng, L. Yue, and Q. Yuan, Adaptive norm selection for regularized image restoration and super-resolution, IEEE Trans. Cyber. 46 (2016) 1388-1399. [17] C. Miao and H. Y. Yu, A General Thresholding Solution for lp (0
45

ACCEPTED MANUSCRIPT [18] T. KÖhler, X. Huang, F. Schebesch, A. Aichert, A. Maier, and J. Hornegger, Robust Multiframe super-resolution employing iteratively Re-weighted minimization, IEEE Trans. Comp. Imaging, 2 (2016) 42-58. [19] B. D. O. Anderson, J. B. Moore, Optimal Filtering, Chap. 3, Prentice-Hall, Inc. Cliffs, New Jersey, 2005. [20] F. Dellaert, C. Tborpe, and S. Thrun, Jacobian image of super-resolved texture maps for model-

Vision, New Jersey,(1998)2-7.

CR IP T

based motion estimation and tracking, in Fourth IEEE Workshop on Application of Computer

[21] S. Farsiu, M. Elad, and P. Milanfar, Video-to-video dynamic superresolution for grayscale and color sequences, EURASIP J. Appl. Signal Process. (2006) 1-15.

[22] C. Newland, D. Gray, and D. Gibbins, Modified Kalman filtering for image super-resolution, in

AN US

Proc. Of the Conf. (IVCNZ06), New Zealand (2006) 79-84.

[23] Z. Wei, F. Tao, and W. Jun, Kalman filter-based method for image superresolution using a sequence of low-resolution images, J. Elec. Imaging, 23(1) (2014) 1-12. [24] V.V. Petrov, Limit Theorems of Probability Theory: Sequences of Independent Random Variables, Clarendon Press, 1995.

M

[25] Z. Wang, Bovik, A.C. Sheikh, E. P., Simoncelli, Image quality assessment: from error visibility to structural similarity, IEEE Trans. on Image Process. 13( 2004) 600-612.

ED

[26] L. Zhang, L. Zhang, X. Mou, and D. Zhang, FSIM: a feature similarity index for image qualtiy assessment, IEEE Trans. on Image Process. 20 (2011) 2378-2386. of

California,

PT

[27]University

Santa

Cruz,

Milanfar

Software

[Online]

Available:

http://users.soe.ucsc.edu/-milanfar/software/sr-datasets.html

CE

[28] S. Zhu , K. Ma, A new diamond search algorithm for fast block matching motion estimation,

AC

IEEE Trans. Image Process. 9 (2000) 287–290.

46