A fractional order variational model for the robust estimation of optical flow from image sequences

Optik 127 (2016) 8710–8727

Contents lists available at ScienceDirect

Optik journal homepage: www.elsevier.de/ijleo

Original research article

A fractional order variational model for the robust estimation of optical flow from image sequences Pushpendra Kumar a,∗ , Sanjeev Kumar a , Balasubramanian Raman b a b

Dept. of Mathematics, Indian Institute of Technology Roorkee, Roorkee 247667, India Dept. of Computer Science & Engineering, Indian Institute of Technology Roorkee, Roorkee 247667, India

a r t i c l e

i n f o

Article history: Received 5 August 2015 Accepted 26 May 2016 Keywords: Fractional derivative Image sequence Optical flow Partial differential equations Variational methods

a b s t r a c t In this paper, a fractional order variational model for estimating the optical flow is presented. In particular, the proposed model generalizes the integer order variational optical flow models. The fractional order derivative describes discontinuous information about texture and edges, and therefore a more suitable in estimating the optical flow. The proposed variational functional is a combination of a global model of Horn and Schunck and the classical model of Nagel and Enkelmann. This formulation yields a dense flow and preserves discontinuities in the flow field and also provides a significant robustness against outliers. The Grünwald–Letnikov derivative is used for solving complex fractional order partial differential equations. The corresponding linear system of equations is solved by an efficient numerical scheme. A detailed stability and convergence analysis is given in order to show the mathematical applicability of the numerical algorithm. Experimental results on various datasets verify the validity of the proposed model. © 2016 Elsevier GmbH. All rights reserved.

1. Introduction Recovery of optical flow from a sequence of images is a key problem in computer vision and image processing. Optical flow is defined as a 2D velocity vector which arises either due to the motion of the objects in the scene or by the motion of the observer/camera. The objective of the problem is to determine the displacement in terms of optical flow between two image frames. Therefore, it is known as correspondence problem in vision system. Optical flow furnishes the dynamic information of an object between the image frames of a scene, i.e., how many units the pixel/object has moved compared to the previous frame. It has a much higher level applications in vision system such as 3D reconstruction, automatic navigation, video surveillance, human action understating and medical diagnosis [1–4]. Extraction of the optical flow has been also used as a preprocessing step in many vision algorithms such as visual control, motion parameter estimation and image segmentation [5–8]. Estimation of the flow field is considered as an ill-posed problem. Thus, it is required to consider some additional constraints in order to regularize the flow field during optical flow estimation. Many models have been proposed to determine the optical flow in the literature starting from the seminal work [9,10] and achieved an impressive level of accuracy. A global approach proposed by Horn and Schunck [9] yields a dense flow, but it is experimentally sensitive to noise. A local approach proposed by Lucas and Kanade [10] is robust under noise, but it is unable to yield dense flow.

∗ Corresponding author. E-mail address: [email protected] (P. Kumar). http://dx.doi.org/10.1016/j.ijleo.2016.05.118 0030-4026/© 2016 Elsevier GmbH. All rights reserved.

P. Kumar et al. / Optik 127 (2016) 8710–8727

8711

In the recent past, differential variational models are one of the most successful approaches for the estimation of optical flow and therefore achieved a good attention of researchers. The reason behind it is their simplicity and advantages in modeling, estimation and quality [11–15]. In order to improve the estimation accuracy of the optical flow models, many assumptions have been made in the variational models to get a significantly better performance. Recent models like [11,12,15–18] include the additional constraints such as gradient constancy assumption and the convex robust data term to obtain the impressive level of performance under outliers, noise and to avoid the problem of local minima. Many of the variational models such as [11–13,19–23] considered different constraints or regularizers to get piecewise optical flow and preserve discontinuity at the object boundaries. In [12,17], the variational functionals were proposed to estimate the optical flow for large motion. Moreover, some researchers have proposed the motion segmentation or parametric model to break the motion field into several smooth and piecewise parts [11,20,24,25]. Nevertheless, these models still have lacked to yield dense and discontinuity preserving optical flow together due to the inherent characteristics of the parametric models. All these models are based on integer order differentiation. These integer order differential variational models can be generalized from integer order to fractional order and categories into the class of fractional order variational models. The fractional order differential based methods are now more popular in image processing applications such as texture analysis, image de-blurring and image restoration [26–34]. The core idea of fractional order differentiation was introduced at the end of sixteen century and later published in the nineteenth [35,36]. The fractional order differential models deal with differentiation of arbitrary order and therefore known by the generalization of integer order variational models. In the literature, numerous authors have been proposed the definitions of fractional order differentiation. But, the quite popular fractional order differentiations among them are Caputo definition [37], Grünwald–Letnikov definition [38] and Riemann–Liouville definition [39]. The salient difference between the fractional and integer order differentiations is that we can find the fractional derivatives even if the function is not continuous whereas integer order derivative failed. Thus, fractional order derivatives efficiently provides the discontinuous information about texture and edges in the optical flow field. Due to non-local property, fractional order derivatives give an excellent tool for the description of memory and modeling of many characteristic phenomena in engineering and sciences. Moreover, fractional order differentiations provide their optical fractional order corresponding to a stable solution. 1.1. Contribution In this paper, we proposed a fractional order variational model to estimate the optical flow from a sequence of images. The novelty of the proposed model compared to previous approaches is in three-fold. First, the proposed model generalizes the existing variational models from integer order to fractional order. Second, the given variational functional is formed from a global model of Horn and Schunck [9] and the classical model of Nagel and Enkelmann [13]. Thus, it combined the advantages of each of them. This helps to lead a dense flow over a region and preserves discontinuities in the optical flow. But, the recent model [28] is based on [9], which provides only dense flow and unable to preserve discontinuities. Third, due to the presence of fractional order derivatives, the proposed model efficiently handles texture and edges. The related Euler–Lagrange equations are obtained corresponding to the fractional order variational model. The numerical discretization of the complex fractional order partial differential equations is based on Grünwald–Letnikov derivative and solved by an efficient numerical scheme. The performance of the proposed model has been tested on various spectrum databases, and compared with some recently published techniques such as [28,40]. The convergence analysis of the used numerical scheme is provided to support the applicability of the algorithm. The robustness of the proposed model are shown under a considerable amount of noise. The rest of the paper is organized in the following sections: In Section 2 some preliminaries related to the fractional order differentiations have been given. Section 3 describes the optical flow constraints followed by proposed fractional variational model for optical flow estimation. Section 4 describes the numerical scheme with its convergence analysis and a pseudo code of the algorithm. Section 5 describes the experimental datasets and evaluation methods followed by the experimental results. Conclusions are given in Section 6. 2. Mathematical preliminaries of fractional order derivatives 2.1. Riemann–Liouville definition Let g(t) ∈ L1 ([a, b]) be a function and ˛ be a fractional value between m − 1 to m, m ∈ Z+ . The left Riemann–Liouville fractional integral of order ˛ is defined as [39] ˛ a It g(t)

=

1  (˛)



t

(t − )˛−1 g()d,

t ∈ [a, b]

(1)

a

where (˛) represents the gamma function, which can be defined as an improper integral

 (˛) = 0

∞

e−t t ˛−1 dt

(2)

8712

P. Kumar et al. / Optik 127 (2016) 8710–8727

Similarly, the right Riemann–Liouville fractional integral of order ˛ is defined as ˛ t Ib g(t) =

1  (˛)



b

( − t)˛−1 g()d,

t ∈ [a, b]

(3)

t

The Riemann–Liouville fractional derivatives of order ˛ of a continuous function g(t) ∈ [a, b] is defined as ˛ a Dt g(t)

= =

dm dt

m−˛ g(t)) m (a It

dm 1  (m − ˛) dt m



(4)

t

(t − )m−˛−1 g()d,

t ∈ [a, b]

a

and ˛ t Db g(t)

 d m

=

−

dt

m−˛ g(t) t Ib



d 1 −  (m − ˛) dt

=

m 

(5)

b

( − t)m−˛−1 g()d,

t ∈ [a, b]

t

The derivatives defined in (4) and (5) are called the left and right Riemann–Liouville fractional derivatives of order ˛, respectively. The Riemann–Liouville definition of fractional derivatives and integrals of order ˛ satisfy the linear and additive property. For ˛ ∈ Z+ , the Riemann–Liouville derivatives provide the similar results as integer order differentiation. 2.2. Grünwald–Letnikov definition The Grünwald–Letnikov derivative [38] of order ˛ is defined as the limit of a sum given by 1  (−1)r h→0 h˛ n

Dt˛ g(t) = lim

r=0

 

  ˛ r

g(t − rh)

(6)

˛ is the binomial coefficient and h is mesh grid size. This derivative formula exhaustively used for discretization r of integer order derivatives.

where

3. Mathematical formulation of variational optical flow model 3.1. Optical flow as variational problem Let I(x, y, t) and I(x + ıx, y + ıy, t + ıt) be the intensities of an image pixel in two consecutive frames having spatiotemporal coordinates (x, y, t) and (x + ıx, y + ıy, t + ıt), respectively. Then by the assumption of [9], I(x + ıx, y + ıy, t + ıt) = I(x, y, t)

(7)

where I :  ⊂ R3 → R represents a rectangular images. By applying Taylor series expansion upto first order in the left-hand side of (7) and simplifying, we get uIx + vIy + It = 0

(8)

The terms Ix , Iy and It denote the partial derivatives of the image function with respect to x, y and t, respectively. The notation u(x, y) and v(x, y) are called the optical flow components in x and y-directions, respectively. The expression derived in equation (8) is called the optical flow constraint (OFC) and plays an important role in optical flow estimation. The solution of optical flow estimation problem is under-determined and require one more constraint on (u, v) in order to determine a unique solution. This constraint is known as smoothness constraint. It tells that neighboring pixels should have the same flow field in a small region of the image. The seminal work of Horn and Schunck [9] minimizes the following energy functional to estimate the optical flow







(uIx + vIy + It )2 +  (|∇ u|2 + |∇ v|2 ) dX

E(u) = 

where  > 0 is a smoothness parameter, X=(x,y) and u = (u, v).

(9)

P. Kumar et al. / Optik 127 (2016) 8710–8727

8713

3.2. Proposed variational optical flow model In order to improve the robustness of the variational model and find a more accurate flow field, a modified version of the energy functional (9) is defined as,







(uIx + vIy + It )2 + ˇ2 (u2 + v2 ) +  (|∇ u|2 + |∇ v|2 ) dX

E(u) =

(10)



where ˇ is a small positive constant. The motivation behind the energy functional (10) is that it combined the global model of Horn and Schunck [9] and classical model of Nagel and Enkelmann [13]. The model of [9] provides dense and smooth flow inside each motion field, but unable to preserve discontinuity in the optical flow and therefore flow propagates in all the directions. The model [13] provides more sharp edges and boundaries in the flow fields and also preserve discontinuity in the optical flow. Thus, it fuses the advantages of each of them. The resulting variational functional (10) provides more accurate dense optical flow and preserves discontinuities as well. 3.3. Proposed fractional order variational optical flow model The fractional order variational model of the above energy functional (10) is defined as







(uIx + vIy + It )2 + ˇ2 (u2 + v2 ) +  (|D˛ u|2 + |D˛ v|2 ) dX

E(u) =

(11)



where D˛ := (Dx˛ , Dy˛ )T represents the left fractional derivative operator of Riemann–Liouville and |D˛ u| =



(Dx˛ u)2 + (Dy˛ u)2 .

Since, fractional order differentiation is a generalization of integer order differentiation, which computes the derivatives of fractional order ˛ ∈ R+ . Therefore, the proposed fractional order variational model can be summarized as: • When ˛ = 1 and ˇ = 0, then the proposed model can be seen as the integer order variational optical flow model proposed by Horn et al. [9]. • When ˛ = 2 and ˇ = 0, then the proposed fractional order optical flow model can be seen as the second order variational optical flow model presented by Yuan et al. [41]. • When ˛ = 1, then the proposed model can be seen as the integer order variational optical flow model proposed by Kumar et al. [42]. • When ˇ = 0, then the proposed model represents the fractional order variational optical flow model given by Chen et al. [28]. Therefore, it is obvious that when ˛ ∈ (0, 1) or (1, 2), the proposed fractional order variational optical flow model can be considered as the generalization of the integer order variational optical flow models. Now, this variational functional (11) has to be minimized in order to estimate the flow vector u = (u, v). Therefore, the above variational functional (11) precisely can be written as







(uIx + vIy + It )2 + ˇ2 (u2 + v2 ) +  (|Dx˛ u|2 + |Dy˛ u|2 + |Dx˛ v|2 + |Dy˛ v|2 ) dX

E(u) =

(12)



In order to find the corresponding Euler–Lagrange equations of (12), let us consider u* (x, y) and v∗ (x, y) are the functions that have to be desired such that for any test functions (x, y) and ϕ(x, y) ∈ C∞ , u and v are defined as [28] u(x, y) = u∗ (x, y) + (x, y) v(x, y) = v∗ (x, y) + ϕ(x, y)

(13)

where  ∈ R. Thus, the Riemann–Liouville fractional derivative of (13) with respect to x and y are Dx˛ u(x, y) = Dx˛ u∗ (x, y) + Dx˛ (x, y) Dy˛ u(x, y) = Dy˛ u∗ (x, y) + Dy˛ (x, y) Dx˛ v(x, y) = Dx˛ v∗ (x, y) + Dx˛ ϕ(x, y) Dy˛ v(x, y) = Dy˛ v∗ (x, y) + Dy˛ ϕ(x, y)

(14)

8714

P. Kumar et al. / Optik 127 (2016) 8710–8727

using Eqs. (13) and (14) in (12), we obtain that the functional (12) is a function of  only for each (x, y) and ϕ(x, y), which is given as





(Ix (u∗ + ) + Iy (v∗ + ϕ) + It )2 + ˇ2 ((u∗ + )2 + (v∗ + ϕ)2 )

E() = 



+  (|Dx˛ u∗ + Dx˛ |2 + |Dy˛ u∗ + Dy˛ |2 + |Dx˛ v∗ + Dx˛ ϕ|2 + |Dy˛ v∗ + Dy˛ ϕ|2 ) dX

(15)

for getting the extremum of the above variational functional (15), differentiate it with respect to  and put  = 0 [43], we get





{(Ix u∗ + Iy v∗ + It ) + (Ix + ϕIy )}(Ix + ϕIy ) + ˇ2 ((u∗ + ) + (v∗ + ϕ)ϕ)

E() = 



+

(Dx˛ u∗ + Dx˛ )Dx˛  + (Dy˛ u∗ + Dy˛ )Dy˛  + (Dx˛ v∗ + Dx˛ ϕ)Dx˛ ϕ + (Dy˛ v∗ + Dy˛ ϕ)Dy˛ ϕ



dX

(16)

and





((Ix u∗ + Iy v∗ + It )(Ix + ϕIy )) + ˇ2 (u∗ + ϕv∗ )

E(0) = 



+

(Dx˛ u∗ )Dx˛  + (Dy˛ u∗ )Dy˛  + (Dx˛ v∗ )Dx˛ ϕ + (Dy˛ v∗ )Dy˛ ϕ



dX

(17)

rewriting equation (17)









{MI x + ˇ2 u∗ +  (Dx˛∗ Dx˛ u∗ + Dy˛∗ Dy˛ u∗ )} + ϕ MI y + ˇ2 v∗ +  (Dx˛∗ Dx˛ v∗ + Dy˛∗ Dy˛ v∗ )

E(0) =

dX

(18)



where D˛* is the right fractional order derivative operator of Riemann–Liouville and M = Ix u∗ + Iy v∗ + It . After putting E (0) = 0 and equating the coefficients of  and ϕ, we get (Ix u∗ + Iy v∗ + It )Ix + ˇ2 u∗ +  (Dx˛∗ Dx˛ u∗ + Dy˛∗ Dy˛ u∗ ) = 0 (Ix u∗ + Iy v∗ + It )Iy + ˇ2 v∗ +  (Dx˛∗ Dx˛ v∗ + Dy˛∗ Dy˛ v∗ ) = 0

(19)

This derivation composes a system of equations for each pixels of the image. 4. Numerical solution 4.1. Derivative discretization scheme Let dimension of our computed flow fields u and v are same as the image sequence of m × n pixels. Then u and v can be discretized as u(i, j) = u(i h, j h),

i = 0, 1, . . ., m − 1

v(i, j) = v(i h, j h),

j = 0, 1, . . ., n − 1

where (i, j) represents the pixel position in an image and h is mess grid size. In image processing, h is chosen to be 1. Now, using the Grünwald–Letnikov definition, the left fractional derivative of order ˛ of u can be discretized as Dx˛ u(i, j) =

∞ 

w(˛)p u(i − p, j)

(20)

p=0

 

where w(˛)p

= (−1)

p

˛ p

, which can be obtained from the coefficients of power series (1 − z)˛ using the following recursive

formula: w(˛)0 = 1,

w(˛)p =



1−

˛+1 p



w(˛)p−1 ,

p = 1, 2, . . .

(21)

using (5) and (20), we get Dx˛∗ Dx˛ u(i, j) =

0 

(˛)

w|p| u(i − p, j) +

p=−∞

∞  p=0

(˛)

wp u(i − p, j)

(22)

P. Kumar et al. / Optik 127 (2016) 8710–8727

Since,

∞

(˛)

p=0

wp

Dx˛∗ Dx˛ u(i, j) =

8715

= 0, therefore (22) can be rewritten as −1 

(˛)

w|p| ∇ u(i − p, j) +

p=−∞

∞ 

(˛)

wp ∇ u(i − p, j)

(23)

p=1

where ∇ u(i − p, j) = u(i − p, j) − u(i, j). For applications point of view in image processing, the above formula (23) can be approximated as Dx˛∗ Dx˛ u(i, j) ≈

−1 

(˛)

w|p| ∇ u(i − p, j) +

p=−W

W 

(˛)

wp ∇ u(i − p, j)

(24)

p=1

where W represents the window mask size. Similarly, derivative of the flow field u with respect to y, we get Dy˛∗ Dy˛ u(i, j) ≈

−1 

(˛)

w|p| ∇ u(i, j − p) +

p=−W

W 

(˛)

wp ∇ u(i, j − p)

(25)

p=1

Combining (24) and (25), we obtain the discrete fractional order derivative formula for u, which is given as



Dx˛∗ Dx˛ u(i, j) + Dy˛∗ Dy˛ u(i, j) ≈

(˛) wp ¯ (u(¯i, ¯j) − u(i, j)) ij

(26)

(¯i,¯j) ∈ (i,j)

where (i, j) is the set of neighbor pixels of pixel position (i, j) in the directions of x and y. Here, pij¯ = max(|¯i − i|, |¯j − j|). In the similar way, we obtain the discrete fractional order derivative formula for v, which is defined as



Dx˛∗ Dx˛ v(i, j) + Dy˛∗ Dy˛ v(i, j) ≈

(˛)

wp ¯ (v(¯i, ¯j) − v(i, j)) ij

(27)

(¯i,¯j) ∈ (i,j)

4.2. Numerical scheme For better understanding, we consider the following abbreviations as Ixx = Ix2 (i, j), Iyy = Iy2 (i, j), Ixy = Ix (i, j)Iy (i, j), Ixt = Ix (i, j)It (i, j), Iyt = Iy (i, j)It (i, j), u¯ = u(¯i, ¯j),

v¯ = v(¯i, ¯j) Using (26) and (27) in (19), we get (Ixx + ˇ2 )u + Ixy v + 



(˛)

wp ¯ (u¯ − u) = −Ixt ij

(¯i,¯j) ∈ (i,j)

Ixy u + (Iyy + ˇ2 )v + 



(˛)

wp ¯ (¯v − v) = −Iyt

(28)

ij

(¯i,¯j) ∈ (i,j)

This system of linear equations (28) is sparse and very large because the number of rows and columns are equal to double of image pixels. We used Gauss–Seidel method [44,45] as an iterative scheme to solve this system (28). The iteration equations for u and v are given as un+1 = u¯ n +

(Ixx − ˇ2 Ixx − ˇ2 Iyy + ˇ2  − ˇ4 )u¯ n + Ixy v¯ n + ( − ˇ2 )Ixt D

(29)

vn+1 = v¯ n +

Ixy u¯ + (−ˇ2 Ixx + ( − ˇ2 )Iyy + ˇ2  − ˇ4 )¯v − (ˇ2 − )Iyt D

(30)

where  =

(˛) w (¯i,¯j) ∈ (i,j) pij¯

and D is the determinant of the system of equations (28).

8716

P. Kumar et al. / Optik 127 (2016) 8710–8727

4.3. Convergence analysis In order to analyze the convergence of the numerical algorithm, we require to determine under what conditions it converges to the exact solution? i.e., under what conditions it produces a sequence of approximations? Let our system of equations in (28) be written in the following form Ax = b



where A =

Ixx + ˇ2 − 

(˛) w (¯i,¯j) ∈ (i,j) pij¯

Ixy Iyy + ˇ2 − 

Ixy

⎡



−

⎢ ⎣ −

b=⎢

(˛)

wp ¯ u¯ − Ixt ij

(¯i,¯j) ∈ (i,j)



(˛) wp ¯ v¯ − Iyt ij

(31)



(˛) w (¯i,¯j) ∈ (i,j) pij¯

  , x=

u

v

and

⎤ ⎥ ⎥ ⎦

(¯i,¯j) ∈ (i,j)

In general, iteration method can be demonstrated in the following form in order to explore the convergence property of the solution of (31) xn+1 = Hxn + B,

n = 0, 1, 2, . . .

(32)

where xn+1 and xn are the approximations of x at the (n + 1)th and nth iterations, respectively. The iteration matrix H =− (D + L)−1 U and vector B = (D + L)−1 b are known. The terms L, U and D are the lower triangular, upper triangular and diagonal matrices decomposed from A such that: A=L+U+D where,

 L=

0

0

Ixy

0



 , U=

0

Ixy

0

0



and

⎡



Ixx + ˇ2 − 

⎢ ⎣0

D=⎢

(˛)

wp ¯

ij

⎤

0

(¯i,¯j) ∈ (i,j)



Iyy + ˇ2 − 

(˛) wp ¯ ij

⎥ ⎥ ⎦

(¯i,¯j) ∈ (i,j)

Now

⎡



Ixx + ˇ2 − 

⎢ ⎣ Ixy

(D + L)−1 = ⎢

(˛)

wp ¯

ij

⎤−1

0

(¯i,¯j) ∈ (i,j)

Iyy

+ ˇ2

−

 (¯i,¯j) ∈ (i,j)

Thus,

 H = −(D + L)−1

0

Ixy

0

0



using the norm of vector in (32) is equivalent to

xn+1 =

Hxn + B

≤ H xn + B

(˛) wp ¯ ij

⎥ ⎥ ⎦

P. Kumar et al. / Optik 127 (2016) 8710–8727

8717

Table 1 Details about the eigenvalues and spectral radius of iteration matrix H. Image sequence

max

(H)

Yosemite Translated tree Rubic cube Heart Face Stillcap Home Mans Grove RubberWhale Urban Venus Army Mequon Teddy Wooden

1.0178 1.0040 1.0294 1.0017 1.0034 1.0218 1.0992 1.0129 1.0070 1.0683 1.0062 1.0137 1.0691 1.0707 1.0528 1.0872

1.0178 1.0040 1.0294 1.0017 1.0034 1.0218 1.0992 1.0129 1.0070 1.0683 1.0062 1.0137 1.0691 1.0707 1.0528 1.0872

max = largest eigenvalue and (H) = spectral radius.

it can be further rewritten as

xn+1

= H xn + B

≤ H ( H

xn−1 + B ) + B

= H 2 xn−1 + H

B + B

= H 2 xn−1 + ( H + 1) B

≤ H 2 ( H

xn−2 + B ) + ( H + 1) B

(33)

= H 3 xn−2 + ( H 2 + H + 1) B

.. . ≤ H n+1 x0 + ( H n + · · · + H + 1) B

If H <1 then (33) reduced to the following form

xn+1 ≤ H n+1 x0 +

1

B

1 − H

(34)

Since, H n + · · · + H +1 is a partial sum of the geometric series, therefore it tend to the following form as ∞ 

H n =

n=0

1 1 − H

(35)

because the term H n+1 → 0 as n → ∞. If H <1, then from Eq. (34), we conclude that H <1 will give convergence. The following statements give the convergence guarantee. For further explanation one can refer [44,45]. Theorem 4.1.

Let H be a square matrix. Then

lim H n = 0

n→∞

if H <1 or iff (H) < 1, where (H) represents the spectral radius of matrix H[44]. Theorem 4.2. A necessary and sufficient condition for the convergence of an iterative method of the form (32) is that the eigenvalues of the iterative matrix satisfy: |i (H)| < 1,

i = 1, . . ., n

The detail information about the eigenvalues and the spectral radius of the iteration matrix for all image sequences are given in Table 1. Table 1 shows that the given algorithm is converging for all image sequences used in the experimental study.

8718

P. Kumar et al. / Optik 127 (2016) 8710–8727

The pseudo code of the proposed algorithm is described in Algorithm 1. Algorithm 1.

Pseudo code of the proposed algorithm for optical flow estimation

Input: I1 , I2 (image frames of a sequence), , ˇ, ˛ and iterations Output: Optical flow vector u=(u,v) 1: Compute the derivatives of I1 , I2 2: In each of the level starts with the coarsest(every variable initialize with the current level in the image) 3: Initialize u with zero at the coarsest level 4: Otherwise update with the current level from the previous level 5: if ˛ = integer value 6: Compute integer order derivatives 7: else 8: Compute fractional order derivatives 9: Compute I2 (X + uo ) and ∇ I2 (X + uo ) using bilinear interpolation 10: for j=1 to iterations 11: Compute u and v using Eqs. (29) and (30) 12: j=j+1 13: end 14: end 15: Return u

5. Experiments, results and discussion 5.1. Datasets In order to assess the performance of any optical flow model, experimental datasets have a major role. In this context, various databases are available online. In this paper, we analyze the performance of the proposed fractional order variational model on different datasets with varying spectrum. All these datasets have different structure and properties such as texture, motion blur, shadow, etc. The details of the datasets which we have used in the experimental study are given as follows: • • • • •

Synthetic dataset: Yosemite with clouds and Translated tree [46]. Real dataset: Heart and Rubic cube [46]. Thermal dataset: Face, Home, Mans and Stillcap [40,47]. Middlebury training dataset: Grove, RubberWhale, Venus and Urban [48]. Middlebury evaluation dataset: Army, Mequon, Teddy and Wooden [48].

All these sample image datasets are shown in Fig. 1. The detailed information about the synthetic and real datasets are available in [46]. The thermal datasets with their brief description are given in [40]. Middlebury training and evaluation datasets are obtained from [48]. 5.2. Evaluation method 5.2.1. Performance measure For assessing the performance of the proposed fractional order variational model, we performed several experiments on the above mentioned datasets. For comparing the proposed model with the other approaches, we applied some presmoothing steps. These pre-smoothing steps are common nowadays to compute optical flow [46]. This pre-smoothing has been done by the convolution of images with a Gaussian kernel of standard deviation n . In order to measure the performance of an optical flow algorithm, many methods have been proposed in the literature. We estimated the angular error(AE) for measuring the performance of the proposed algorithm. This AE is the angle between the correct flow vector (uc , vc , 1) and the estimated flow vector (ue , ve , 1) (see [46] for details) and defined as,



AE = cos

−1



uc ue + vc ve + 1



(u2c + v2c + 1)(u2e + v2e + 1)

(36)

Arbitrary scaling constant 1 in AE is used as a unit converter from pixels to degrees. It plays an important role in case of zero flows to avoid the problem of dividing by zero. 5.2.2. Statistics The validity of the proposed model is also verified using the following statistical measures: • Average angular error (AAE) • Percentage error (PE) • Standard deviation (STD)

P. Kumar et al. / Optik 127 (2016) 8710–8727

8719

Fig. 1. Sample original image datasets. (a) Yosemite with clouds, (b) Translated tree, (c) Rubic cube, (d) Heart, (e) Face, (f) Stillcap, (g) Home, (h) Mans, (i) Venus, (j) Urban, (k) RubberWhale, (l) Grove, (m) Army, (n) Mequon, (o) Teddy and (p) Wooden.

All these terms are briefly defined in [46]. The smaller the values of these terms will show the higher estimation accuracy of the model. 5.3. Experimental discussion All the experiments of the proposed algorithm have been performed on the MATLAB R2011b platform (64 bit version) for windows of 4 GB RAM. We performed six experiments for a detailed evaluation, analysis and comparisons of the proposed algorithm with the existing algorithms. The default values of the smoothing parameter  and ˛ are 100 and 1, respectively. A window mask of size 3 × 3 is used to find the fractional order derivatives in all the experiments. Qualitative and quantitative experimental results are provided to demonstrate the validity of the proposed fractional order variational model instead of integer order variational model. The qualitative results are given in the form of vector representation and color maps of the optical flow field. The quantitative results are shown in terms of statistical results. In the first experiment, we have estimated the statistical results corresponding to fractional order ˛ of the proposed model for all datasets. This relationship between statistical results and ˛ are shown in Figs. 2 and 3. The smaller is the statistical result value provides the better estimation accuracy. Therefore, the optimal value of fractional order ˛ is depend on statistical errors. The obtained results given in Figs. 2 and 3 provide different optimal fractional order for different datasets. The optimal value of fractional order is corresponds to the stable solution. In the second experiment, we have estimated the optical flow results for synthetic and real image sequences. These image sequences contain several image degradation factors such as camera motion, rotation of objects and illumination effect, etc. The qualitative results are shown in Fig. 4. The quantitative (statistical) results for Yosemite with clouds are

8720

P. Kumar et al. / Optik 127 (2016) 8710–8727

Fig. 2. Optical flow quantitative results (AAE and STD) of the proposed model corresponding to fractional order ˛ for given image sequences: (a) Yosemite with clouds, (b) Translated tree, (c) Rubic cube, (d) Heart, (e) Face, (f) Stillcap, (g) Home and (h) Mans. Here NAN = Bigger value or Not a number.

compared with the other existing techniques in Table 2. Table 2 validates that the proposed fractional order variational model gives comparatively better results than integer order variational models. Qualitative results show that the model preserves discontinuity due to its fractional order derivatives, which verifies the importance of the proposed model. In the third experiment, we have estimated the optical flow results for the thermal image sequences. Since, this dataset is an important part of the vision system and has many applications in tracking and sports. The importance of these datasets is that their intensity values do not effected by shadow or lights. Qualitative results for these datasets are shown in Fig. 5 and quantitative results are compared with the results of algorithm [40] in Fig. 6. This comparison verifies the significance of the proposed model. In the fourth experiment, we analyzed the performance of the proposed algorithm on Middlebury training image datasets. All these datasets have different structures and contain many image degradation components such as several independent motion, texture, substantial motion discontinuities, motion blur, etc. The obtained qualitative results are shown in Fig. 7 and

P. Kumar et al. / Optik 127 (2016) 8710–8727

8721

Fig. 3. Optical flow quantitative results (AAE and STD) of the proposed model corresponding to fractional order ˛ for given Middlebury image sequences: (a) Grove, (b) RubberWhale, (c) Urban, (d) Venus, (e) Army, (f) Mequon, (g) Teddy and (h) Wooden.

Table 2 Comparisons between the results of proposed model and the results from the literature corresponding to Yosemite image sequence with clouds. Yosemite with clouds Technique

AAE

STD

Anandan [11] Nagel [13] Quenot [49] Sun [50] Fleet [51] Proposed scheme

16.37◦ 12.70◦ 9.93◦ 9.21◦ 5.28◦ 1.3124◦

13.46◦ 16.68◦ 16.16◦ 16.16◦ 14.33◦ 0.1254◦

STD, standard deviation; AAE, average angular error.

8722

P. Kumar et al. / Optik 127 (2016) 8710–8727

Fig. 4. Optical flow results for synthetic and real image datasets. First and second row represent the vector forms and color maps of the estimated optical flows corresponding to the image datsets, respectively.

Fig. 5. Optical flow results for thermal image datasets. First and second row represent the vector forms and color maps of the estimated optical flows corresponding to the thermal image sequences, respectively.

quantitative results are compared with the results of algorithm [9] in Table 3. These results demonstrate that proposed model provides more accurate flow fields and strongly handle the textures and edges. This performance verifies the applicability of the proposed model. In the fifth experiment, we have estimated the optical flow results for Middlebury evaluation image datasets. These image sequences also contain textures, shadow and motion blur. The estimated results are demonstrated in Fig. 8. The quantitative results are compared with results of algorithm [9] in Table 3. These comparisons show the superiority of the

Fig. 6. Comparison of quantitative results of the proposed fractional order model (Proposed) with the results of algorithm (Kumar [40]) for thermal datasets.

P. Kumar et al. / Optik 127 (2016) 8710–8727

8723

Fig. 7. Optical flow results for Middlebury training datasets. First and second row represent the vector forms and color maps of the estimated optical flows corresponding to the image sequences, respectively.

Table 3 Comparisons of quantitative results of the proposed model with the results of algorithm [9] corresponding to given sample image datasets. Algorithms

Horn and Schunck [9]

Image sequence

STD

AAE

Proposed model STD

AAE

Yosemite Translated tree Rubic Heart Face Stillcap Home Mans Grove RubberWhale Urban Venus Army Mequon Teddy Wooden

0.1351 0.0389 0.0802 0.0550 0.2491 0.1224 0.1206 0.0818 0.0630 0.1151 0.0998 0.1564 0.0833 0.9295 0.8258 0.7786

1.3251 1.2774 – – – – – – 1.2325 1.5632 1.1336 1.2322 – –

0.1254 0.0361 0.0737 0.0463 0.0017 0.1144 0.0018 0.0011 0.0561 0.1094 0.0960 0.1466 0.0013 0.0110 0.0097 0.0096

1.3124 1.2678 – – – – – – 1.2188 1.5606 1.1220 1.2211 – –

–

– –

Fig. 8. Optical flow results for Middlebury evaluation datasets. First and second row represent the vector forms and color maps of the estimated optical flows corresponding to the image sequences, respectively.

8724

P. Kumar et al. / Optik 127 (2016) 8710–8727

Fig. 9. Quantitative results (AAE and STD) of the proposed model with the addition of Gaussian noise of standard deviation ( n ) for their optimal fractional order ˛: (a) Synthetic datasets: Yosemite and Translated tree and (b) Real datasets: Rubic cube and Heart.

proposed model. The obtained color maps of the optical flow field are dense and smooth inside a region and also provide discontinuous information about the boundaries. This behavior justifies the accuracy of the proposed model. In the next experiment, we have validated the robustness of the proposed model in the presence of noise. The experiments have been performed by the addition of Gaussian noise with different amount of standard deviation n and zero mean on all the image datasets. The performance has been evaluated in terms of statistical results. The obtained results are shown in Figs. 9–11. These results have found corresponding to their optimal value of fractional order ˛. These figures demonstrate that the performance of the proposed model does not influence more by the addition of Gaussian noise. This behavior justifies the robustness of the proposed model. In addition to it, we compared the quantitative results of the proposed model for Middlebury datasets with the results of the recently published state-of- the-art technique [28] in Fig. 12. We also compared the quantitative results for all datasets with the results of the algorithm [9] in Table 3. A summary of quantitative results without the addition of noise for all datasets corresponding to their optimal fractional order ˛ are given in Table 4. Table 4 Summary of quantitative results of the proposed model for given sample image sequences corresponding to their optimal fractional order ˛. Image sequence

STD

AAE

AE

PE

Yosemite Translated tree Rubic Heart Face Stillcap Home Mans Grove RubberWhale Urban Venus Army Mequon Teddy Wooden

0.1254 0.0361 0.0737 0.0463 0.0017 0.1144 0.0018 0.0011 0.0561 0.1094 0.0960 0.1466 0.0013 0.0110 0.0097 0.0096

1.3124 1.2678 – – – – – – 1.2188 1.5606 1.1220 1.2211 – – – –

0.0644 0.1075 – – – – – – 0.0899 0.0689 0.1744 0.0865 – – – –

0.0372 0.0540 – – – – – – 0.0236 0.0546 0.0241 0.0229 – – – –

AE, average error; PE, percentage error.

P. Kumar et al. / Optik 127 (2016) 8710–8727

8725

Fig. 10. Quantitative results (AAE and STD) of the proposed model with the addition of Gaussian noise of standard deviation ( n ) for Middlebury datasets corresponding to their optimal fractional order ˛: (a) Training datasets and (b) Evaluation datasets.

Fig. 11. Quantitative result (STD) of the proposed model with the addition of Gaussian noise of standard deviation ( n ) for thermal datasets corresponding to their optimal fractional order ˛.

Fig. 12. Comparison of quantitative results of the proposed model (Proposed) with the results of model (Chen [28]) for Middlebury datasets.

8726

P. Kumar et al. / Optik 127 (2016) 8710–8727

All these quantitative and qualitative results show that the proposed model provides comparatively better results than other variational models. The proposed model also provides the discontinuous information about texture and edges at object boundaries whereas integer order fail. The flow obtained for all datasets are smooth and dense over a region and sharp at the boundaries. Thus, this model can be incorporated in the applications of image segmentation as a preprocessing step such as medical imaging and optical character recognition. 6. Conclusions and future work In this paper, a fractional order variational model was proposed for estimating the optical flow. It generalizes the existing variational models from integer to fractional order and provides higher order differential variational models for integer fractional order. The corresponding complex fractional order partial differential equations to the variational model were discretized using the Grünwald–Letnikov derivative and solved by an efficient numerical scheme. The performance of the proposed model analyzed on various datasets of varying spectrum contain texture, shadows and motion blur. The optimal fractional order for which the solution is stable was provided for each image sequence. Experimental results in terms of qualitative and quantitative have been presented. The comparisons of these results with some recently published models show the accuracy and robustness of the proposed model. The convergence analysis of the iterative scheme shows that the algorithm leads to convergence. The proposed model provides dense flow over a region and preserves discontinuity in the flow field, and also efficiently handled texture and edges. This behavior can be justified by the estimated color maps of the optical flow. Future studies involve the proposed fractional order variational model can be extended to the other variational models. The algorithm with its extension can be used to handle color images. Due to discontinuity preserving nature, this model can be incorporated into the applications of image segmentation. Acknowledgement One of the authors, Pushpendra Kumar gratefully acknowledges to Council of Scientific and Industrial Research (CSIR), New Delhi, Government of India for the financial support to carry out this research work. References [1] J. Lan, J. Li, G. Hu, B. Ran, L. Wang, Vehicle speed measurement based on gray constraint optical flow algorithm, Optik – Int. J. Light Electr. Optics 125 (2014) 289–295. [2] H. Kim, K. Sohn, 3D reconstruction from stereo images for interactions between real and virtual objects, Signal Process.: Image Commun. 20 (2005) 61–75. [3] S.J. Lee, G. Shah, A.A. Bhattacharya, Y. Motai, Human tracking with an infrared camera using a curve matching framework, EURASIP J. Adv. Signal Process. 2012 (2012) 1–15. [4] Y. Motai, S.K. Jha, D. Kruse, Human tracking from a mobile agent: optical flow and kalman filter arbitration, Signal Process.: Image Commun. 27 (2012) 83–95. [5] S. Dhou, Y. Motai, Scale-invariant tracking predictor using a pan-tilt-zoom camera, Robotica (2015) 1–25. [6] S.J. Lee, Y. Motai, H. Choi, Tracking human motion with multichannel interacting multiple model, IEEE Trans. Ind. Inform. 9 (2013) 1751–1763. [7] S. Kumar, S. Kumar, N. Sukavanam, B. Raman, Human visual system and segment-based disparity estimation, AEU – Int. J. Electr. Commun. 67 (2013) 372–381. [8] Y. Xu, X. Hu, S. Peng, Blind motion deblurring using optical flow, Optik – Int. J. Light Electr. Opt. 126 (2015) 87–94. [9] B. Horn, B. Schunck, Determining optical flow, Artif. Intell. 17 (1981) 185–203. [10] B. Lucas, T. Kanade, An iterative image registration technique with an application to stereo vision, in: Seventh International Joint Conference on Artificial Intelligence, vol. 81, 1981, pp. 674–679. [11] M.J. Black, P. Anandan, The robust estimation of multiple motions: parametric and piecewise smooth flow, Comput. Vis. Image Underst. 63 (1) (1996) 75–104. [12] T. Brox, A. Bruhn, N. Papenberg, J. Weickert, High accuracy optical flow estimation based on a theory for warping, in: European Conference on Computer Vision (ECCV 2004), vol. 4, 2004, pp. 25–36. [13] H.H. Nagel, W. Enkelmann, An investigation of smoothness constraints for the estimation of displacement vector fields from image sequences, IEEE Trans. Pattern Anal. Mach. Intell. 8 (1986) 565–593. [14] D. Fortun, P. Bouthemy, C. Kervrann, Optical flow modeling and computation: a survey, Comput. Vis. Image Underst. 134 (2015) 1–21. [15] H. Zimmer, A. Bruhn, J. Weickert, L. Valgaerts, A. Salgado, B. Rosenhahn, H.P. Seidel, Complementary optic flow, in: Energy Minimization Methods in Computer Vision and Pattern Recognition, 2009, pp. 207–220. [16] P. Kumar, S. Kumar, B. Raman, A vision based motion estimation in underwater images, in: International Conference on Advances in Computing, Communications and Informatics, 2015, pp. 1179–1184. [17] L. Alvarez, J. Weickert, J. Sánchez, Reliable estimation of dense optical flow fields with large displacements, Int. J. Comput. Vis. 39 (2000) 41–56. [18] J. Weickert, C. Schnörr, A theoretical framework for convex regularizers in PDE-based computation of image motion, Int. J. Comput. Vis. 45 (2001) 245–264. [19] Y. Lei, L. Jinzong, L. Dongdong, Discontinuity-preserving optical flow algorithm, J. Syst. Eng. Electr. 18 (2007) 347–354. [20] E. Mémin, P. Pérez, Hierarchical estimation and segmentation of dense motion fields, Int. J. Comput. Vis. 46 (2002) 129–155. [21] J. Weickert, On discontinuity-preserving optic flow, in: Proceedings of Computer Vision and Mobile Robotics Workshop, 1998. [22] L. Xu, J. Jia, Y. Matsushita, Motion detail preserving optical flow estimation, IEEE Trans. Pattern Anal. Mach. Intell. 34 (2012) 1744–1757. [23] J. Weickert, C. Schnörr, Variational optic flow computation with a spatio-temporal smoothness constraint, J. Math. Imaging Vis. 14 (2001) 245–255. [24] L. Alvarez, R. Deriche, T. Papadopoulo, J. Sánchez, Symmetrical dense optical flow estimation with occlusions detection, Int. J. Comput. Vis. 75 (2007) 371–385. [25] R. Deriche, P. Kornprobst, G. Aubert, Optical-flow estimation while preserving its discontinuities: a variational approach, in: Recent Developments in Computer Vision, 1996, pp. 69–80. [26] K. Kashu, Y. Kameda, M. Narita, A. Imiya, T. Sakai, Continuity order of local displacement in volumetric image sequence, in: Biomedical Image Registration, 2010, pp. 48–59.

P. Kumar et al. / Optik 127 (2016) 8710–8727

8727

[27] D. Chen, Y. Chen, H. Sheng, Fractional variational optical flow model for motion estimation, in: Proceedings of FDA’1, 2010. [28] D. Chen, H. Sheng, Y. Chen, D. Xue, Fractional-order variational optical flow model for motion estimation, Philos. Trans. R. Soc. London A: Math. Phys. Eng. Sci. 371 (2013) 20120148. [29] J. Bai, X.C. Feng, Fractional-order anisotropic diffusion for image denoising, IEEE Trans. Image Process. 16 (2007) 2492–2502. [30] D. Chen, Y. Chen, D. Xue, Fractional-order total variation image restoration based on primal-dual algorithm, Abstr. Appl. Anal. (2013). [31] D. Chen, Y. Chen, D. Xue, Fractional-order total variation image denoising based on proximity algorithm, Appl. Math. Comput. 257 (2015) 537–545. [32] Y.F. Pu, J.L. Zhou, X. Yuan, Fractional differential mask: a fractional differential-based approach for multiscale texture enhancement, IEEE Trans. Image Process. 19 (2010) 491–511. [33] P. Kumar, S. Kumar, B. Raman, A fractional order variational model for tracking the motion of objects in the applications of video surveillance, in: Eighth International Conference on Advanced Computational Intelligence, 2016, pp. 117–123. [34] S.G. Bardeji, I.N. Figueiredo, E. Sousa, Optical Flow with Fractional Order Regularization: Variational Model and Solution Method, 2015, arXiv:1512.01398. [35] K.B. Oldham, The Fractional Calculus, Elsevier, 1974. [36] K.S. Miller, B. Ross, An Introduction to the Fractional Calculus and Fractional Differential Equations, Wiley, New York, 1993. [37] M. Caputo, Linear models of dissipation whose q is almost frequency independent – ii, Geophys. J. Int. 13 (1967) 529–539. [38] K.S. Miller, Derivatives of noninteger order, Math. Mag. (1995) 183–192. [39] B. Riemann, Versuch einer allgemeinen auffassung der integration und differentiation, Gesammelte Werke 62 (1876). [40] S. Kumar, S. Kumar, B. Raman, A variational approach for optical flow estimation in infra-red or thermal images, in: Second International Conference on Image Information Processing, 2013, pp. 56–61. [41] J. Yuan, C. Schörr, G. Steidl, Simultaneous higher-order optical flow estimation and decomposition, SIAM J. Sci. Comput. 29 (2007) 2283–2304. [42] P. Kumar, S. Kumar, A modified variational functional for estimating dense and discontinuity preserving optical flow in various spectrum, AEU – Int. J. Electr. Commun. 70 (2016) 289–300. [43] T.M. Apostol, Calculus, vol. I, John Wiley & Sons, 2007. [44] M. Jain, S. Iyengar, R. Jain, Numerical Methods for Scientific and Engineering Computation, New Age International (P) Limited, 1985. [45] R.W. Hamming, Numerical Methods for Scientists and Engineers, Dover Publications, 2012. [46] J.L. Barron, D.J. Fleet, S. Beauchemin, Performance of optical flow techniques, Int. J. Comput. Vis. 12 (1994) 43–77. [47] J.W. Davis, V. Sharma, Background-Subtraction using Contour-based Fusion of Thermal and Visible Imagery, Comput. Vis. Image Underst. 106 (2007) 162–182. [48] S. Baker, D. Scharstein, J.P. Lewis, S. Roth, M.J. Black, R. Szeliski, A database and evaluation methodology for optical flow, Int. J. Comput. Vis. 92 (2011) 1–31. [49] G.M. Quenot, Computation of optical flow using dynamic programming, in: IAPR Workshop on Machine Vision Applications, Tokyo, Japan, 1996, pp. 249–252. [50] C. Sun, Fast optical flow using 3D shortest path techniques, Image Vis. Comput. 20 (2002) 981–991. [51] D.J. Fleet, A.D. Jepson, M. Jenkin, Phase-based disparity measurement, CVGIP: Image Underst. 53 (2) (1991) 198–210.