Bayesian Optical Flow with Uncertainty Quantification

Ba y esian Optical Flo w with Uncertain t y Quan tification Jie Sun Departmen t of Mathematics, Clarkson Univ ersity , P otsdam, New Y ork 13699, USA Departmen t of Ph ysics, Clarkson Universit y , P otsdam, New Y ork 13699, USA Departmen t of Computer Science, Clarkson Univ ersity , P otsdam, New Y ork 13699, USA Clarkson Cen ter for Complex Systems Science ( C 3 S 2 ), Clarkson Univ ersity , P otsdam, New Y ork 13699, USA E-mail: sunj@clarkson.edu F ernando J. Quev edo Departmen t of Mec hanical and Aerospace Engineering, Clarkson Universit y , P otsdam, New Y ork 13699, USA Clarkson Cen ter for Complex Systems Science ( C 3 S 2 ), Clarkson Univ ersity , P otsdam, New Y ork 13699, USA Erik Bollt Departmen t of Mathematics, Clarkson Univ ersity , P otsdam, New Y ork 13699, USA Departmen t of Electrical and Computer Engineering, Clarkson Univ ersity , Potsdam, New Y ork 13699, USA Clarkson Cen ter for Complex Systems Science ( C 3 S 2 ), Clarkson Univ ersity , P otsdam, New Y ork 13699, USA Abstract. Optical flo w refers to the visual motion observ ed b et w een t wo consecutiv e images. Since the degree of freedom is t ypically m uc h larger than the constrain ts imp osed b y the image observ ations, the straightforw ard formulation of optical flo w as an inv erse problem is ill-posed. Standard approac hes to determine optical flo w rely on form ulating and solving an optimization problem that contains both a data fidelit y term and a regularization term, the latter effectiv ely resolves the otherwise ill-p osedness of the inv erse problem. In this work, we depart from the deterministic formalism, and instead treat optical flow as a statistical in verse problem. W e discuss ho w a classical optical flo w solution can be interpreted as a p oint estimate in this more general framew ork. The statistical approac h, whose “solution” is a distribution of flo w fields, which we refer to as Ba yesian optical flow, allows not only “p oin t” estimates (e.g., the computation of a v erage flow field), but also statistical estimates (e.g., quan tification of uncertaint y) that are b eyond an y standard metho d for optical flo w. As application, we b enchmark Bay esian optical flow together with uncertaint y quan tification using sev eral types of prescrib ed ground-truth flow fields and images. Bayesian Optic al Flow with Unc ertainty Quantific ation 2 1. In tro duction Optical flow reflects the visual motion b et ween consecutiv e images. Determination of optical flo w is imp ortant for applications ranging from machine learning and computer vision [ 28 ], artificial intelligence and robotics [ 4 , 25 ], to scien tific applications from o ceanograph y to w eather forecasting [ 6 , 5 , 10 , 19 , 20 ], to name a few. In classical approac hes optical flow is determined by solving a v ariational optimization problem whic h requires inclusion of a regularization term in additional to data fitting for the problem to b e w ell-p osed [ 1 , 16 , 35 ]. The c hoice of regularization parameter turns out critical for the satisfactory inference of optical flow. Despite the many (comp eting) metho ds of selecting the regularization parameter, none seems to b e most “natural” comparing to the others [ 15 , 30 , 34 ]. Another imp ortant feature of classical optical flow approac hes is that they pro duce a single solution (b y design) as a “p oint estimate”. In practice, the magnitude of the flo w as well as measurement error and noise can v ary significantly from one part of the images to another, and thus the at a given pixel the inferred flo w v ector can b e asso ciated with either large or small uncertain ty . Although not captured in classical optical flow, suc h “uncertaint y” would provide v aluable information if attainable as part of the solution. In this work w e adopt a statistical in version framework for the computation of optical flow. In this framework the estimation of optical flow is reform ulated from the Ba yesian p ersp ectiv e as a statistical in v ersion problem. F rom this new, statistical p ersp ectiv e, different t yp es of information black and mo del assumptions are collectiv ely fused to give rise to a p osterior distribution of the unkno wns of in terest. The p osterior distribution, whic h is t ypically sampled using some appropriately designed Marko v c hain Mon te Carlo sc heme, can b e further used to deriv e v arious estimates. Imp ortantly , unlik e the point estimate of optical flo w obtained b y classical v ariational approaches, the prop osed statistical inv ersion approac h is a metho dological w ay of uncertaint y propagation to pro duce a distribution of candidate optical flo w fields from which v arious statistical properties can be extracted, including ensemble av erage estimation and uncertain ty quantification. The rest of the pap er is organized as follows. In Section 2 w e review the standard optical flo w setup and discuss how it can b e treated as an inv erse problem defined in finite dimensions. In Section 3 we w e review the basics of in verse problems, including the standard least squares approac h and Tikhono v regularization. Then, in Section 4 w e presen t a statistical inv ersion approach for optical flow, and discuss the choice of priors, mo dels, sampling pro cedure and computational algorithms. In Section 5 w e sho w case the utilit y statistical in version approac h for optical flo w using sev eral b enc hmark examples of flo w fields and noisy image data. Finally , Section 6 contains a conclusion and discussion of sev eral unsolv ed and unattempted issues that might b e of in terest for future research. Bayesian Optic al Flow with Unc ertainty Quantific ation 3 2. Optical Flow as an In verse Problem Giv en a sequence of consecutively captured images, the visual relative motion b etw een them is commonly referred to as optical flo w, which can often provide insigh t ab out the actual ph ysical motion. The inference of optical flo w is an outstanding scien tific question, whic h requires making assumptions about the underlying motion as well as the measuremen t pro cess. 2.1. Pr oblem Setup Consider tw o single-channeled (typically gra yscale) digital images (“pictures”) taken from the same scene at tw o nearb y time instances. The image data are represented by t wo matrices F = [ F ij ] n x × n y and G = [ G ij ] n x × n y . Th us each image contains n x × n y pixels, defined on a common tw o-dimensional subspace Ω. The goal of optical flow estimation is to infer a flo w field, defined b y matrices U = [ U ij ] n x × n y and V = [ V ij ] n x × n y where h U ij , V ij i represents the optical “velocity” o ccurring at the ( i, j )-th pixel inferred from the tw o images. The image data F and G are often regarded as sampled data from smo oth functions F ( x, y ) and G ( x, y ), with F ij = F ( x i , y j ) and G ij = G ( x i , y j ) where { ( x i , y j ) } ( i =1 ,...,n x ; j =1 ,...,n y ) are grid p oin ts from a spatial domain Ω. Thus U and V can be view ed as discrete spatial samples of a smo oth 2D flo w field − → W ( x, y ) = h U ( x, y ) , V ( x, y ) i , whic h is defined on Ω that captures the visual optical motion o ccurring b etw een the tw o observ ed images. 2.2. V ariational Appr o ach of Inferring Optic al Flow The classical v ariational approac h of optical flow starts by defining an “energy” functional whose minimization yields an estimation of the optical flo w field [ 1 ]. One of the most widely used functional was prop osed b y Horn and Sc hunc k in 1981 [ 16 ], giv en by E ( U, V ) = Z Z Ω ( F x U + F y V + F t ) 2 dxdy + α Z Z Ω ( k∇ U k 2 + k∇ V k 2 ) dxdy , (1) where U ( x, y ) and V ( x, y ) are smo oth functions defined o ver Ω whic h represen t a candidate flo w field. In the Horn-Sch unck functional, the first term is often referred to as data fidelity as it measures the deviation of the total image intensit y from b eing conserv ative, that is, ho w muc h do es the mo del fit deviates from the condition dF /dt = F x U + F y V + F t = 0 . (2) The second term measures the solution r e gularity by p enalizing solutions that hav e large spatial gradien ts and is called the regularization term. The relativ e emphasis of smo othness as compared to “fitting” the total image intensit y conserv ation equation ( 2 ) is con trolled b y the positive scalar α whic h is called a r e gularization p ar ameter . The main role of the regularization term is to ensure that the minimization of the functional is a wel l p osed problem. Without the regularization term, the problem is ill-p osed. Bayesian Optic al Flow with Unc ertainty Quantific ation 4 Giv en α , the functions U and V that minimize the Horn-Sc h unck functional ( 1 ) satisfy the Euler-Lagrange equations ( F x ( F x U + F y V + F t ) = α ( U xx + U y y ) , F y ( F x U + F y V + F t ) = α ( V xx + V y y ) , (3) whic h are typically solv ed by some iterativ e sc heme ov er a finite set of spatial grid p oin ts [ 16 , 29 ] to pro duce an estimation of the optical flo w. Alternativ ely , one could also discretize the functional ( 1 ) itself to yield a finite-dimensional inv erse problem as discussed in Section 2.3 b elow with solution strategy reviewed in Section 3 . 2.3. Finite-dimensional r epr esentation of the V ariational Optic al Flow F unctional As discussed in the previous section, the classical v ariational approach of optical flo w w orks by first form ulating and minimizing a functional o ver smo oth v ector fields, and then ev aluating the obtained v ector field at the grid p oints on whic h the original image data are giv en. Here w e approach the problem from a differen t route, by first discretizing the functional ( 1 ) to conv ert the functional minimization (an infinite- dimensional problem) in to a finite-dimensional linear inv erse problem, and then solving the in v erse problem to yield solutions which giv e the v alues of a vector field defined ov er the same grid p oin ts as the image data. The remainder of this section will b e fo cused on the con v ersion of the functional ( 1 ) in to a finite-dimensional function defined on a set of uniformly distributed grid p oin ts { ( x i , y j ) } i =1 , 2 ,...,n x ; j =1 , 2 ,...,n y , (4) where x i +1 − x i = ∆ x and y j +1 − y j = ∆ y are the spacing in the x direction and the y direction, resp ectively . The con v ersion will b e ac hieved b y approximating the in tegrals in Eq. ( 1 ) with appropriately derived summations ov er the grid p oin ts. F or notational con venience, we use a b old lo wercase v ariable to denote the vectorization of a matrix. F or example, the b oldface v ector q denotes the column vector obtained b y “v ertically stac king” the columns of a matrix Q = [ − → Q 1 , − → Q 2 , . . . , − → Q n ] in order [ 12 ], where − → Q i denotes the i -th column of Q . That is, q = vec( Q ) =      − → Q 1 − → Q 2 . . . − → Q n      . (5) First let us consider the data fidelity term: R R Ω ( F x U + F y V + F t ) 2 dxdy . The spatial deriv atives F x ( x, y ) and F y ( x, y ) can b e approximated b y a finite difference scheme. F or example, the simple forwar d differ enc e yields the appro ximations: ( F x ( x, y ) ≈ 1 ∆ x [ F ( x + ∆ x, y ) − F ( x, y )] , F y ( x, y ) ≈ 1 ∆ y [ F ( x, y + ∆ y ) − F ( x, y )] . (6) Bayesian Optic al Flow with Unc ertainty Quantific ation 5 W e next express these deriv ativ es as op erations on the column v ector f . T o do this, we define matrix S k = [ S ( ij ) k ] k × k as S ( ij ) k = ( − δ ij + δ i +1 ,j , if i < k ; − δ i,j − 1 + δ i,j , if i = k . (7) The forw ard difference applied to f can b e represen ted as ( f x ≈ Q x f , f y ≈ Q y f , (8) where ( Q x ≡ 1 ∆ x [ I n ⊗ S m ] , Q y ≡ 1 ∆ y [ S n ⊗ I m ] . (9) The temp oral deriv ativ e can b e estimated from the data b y the direct pixel-wise difference b et ween the tw o images, to yield f t ≈ g − f . (10) With these definitions and approximations, w e obtain a discretized version of the conserv ation equation ( 2 ) expressed as a finite-dimensional linear system: A x = b , (11) where      A = [diag( f x ) , diag( f y )] , x = [ u > , v > ] > , b = − f t . (12) Here diag( f x ) and diag( f y ) represent the diagonal matrices whose diagonal elemen ts are given by the entries of the vector f x and f y , resp ectively . F rom this connection we appro ximate the first integral in the functional ( 1 ) as k A x − b k 2 where k · k denotes the standard Euclidean norm. Next w e develop a finite-dimensional approximation of the regularization term in the functional ( 1 ). This requries discretization of ∇ U and ∇ V . Using a similar forward difference to appro ximate the partial deriv ativ es, w e obtain ( ∇ U ≈ 1 ∆ x [ U ( x + ∆ x, y ) − U ( x, y )] + 1 ∆ y [ U ( x, y + ∆ y ) − U ( x, y )] , ∇ V ≈ 1 ∆ x [ V ( x + ∆ x, y ) − V ( x, y )] + 1 ∆ y [ V ( x, y + ∆ y ) − V ( x, y )] . (13) F or the vectorized v ariables u and v , we hav e ( ∇ u = u x + u y ≈ ( Q x + Q y ) u , ∇ v = v x + v y ≈ ( Q x + Q y ) v , (14) where Q x and Q y are defined in Eq. ( 9 ). Consequen tly , we obtain the approximation of the second in tegral in the Horn-Sch unc k functional ( 1 ) as ( R R Ω k∇ U k 2 dxdy ≈ u > x u x + u > y u y ≈ u > [ Q > x Q x + Q > y Q y ] u , R R Ω k∇ V k 2 dxdy ≈ v > x u y + v > y v y ≈ v > [ Q > x Q x + Q > y Q y ] v , (15) Bayesian Optic al Flow with Unc ertainty Quantific ation 6 whic h then giv es Z Z Ω ( k∇ u k 2 + k∇ v k 2 ) dxdy ≈ x > Q x , (16) where the matrix Q = I 2 ⊗ [ Q > x Q x + Q > y Q y ] . (17) Therefore, the v ariational optical flo w formulation ( 1 ) can b e reformulated at a finite spatial resolution as an inv erse problem, where, for a given parameter v alue α , the corresp onding solution is giv en b y solving the following (regularized) optimization: min x  k A x − b k 2 + α x > Q x  . (18) This is a standard approach in inv erse problems, form ulated as a least squares problem with Tikhono v regularization, with more details to b e presented in the next section. 3. Inv erse Problem in Finite Dimensions Although there has b een a great deal of progress on the mathematical characterization of in verse problems in the field of functional analysis, a practical problem often concerns finding a solution in a finite-dimensional space. A t a fundamental level, the most common in verse problem stems from a linear mo del [ 15 , 17 , 34 ] b = A x + η , (19) where b ∈ R m is a column v ector of observe d data , A = [ a ij ] m × n ∈ R m × n is a (known) matrix representing the underlying mo del, the column v ector η ∈ R n denotes (additive) noise, and x ∈ R n is the vec tor of unknowns to b e inferred. Giv en A and b , the problem of inferring x in Eq. ( 19 ) is called an inverse pr oblem b ecause rather than direct “forw ard” computation from the mo del, it requires a set of indirect, “backw ard”, or “in v erse” op erations to determine the unkno wns [ 34 ]. Dep ending on the rank and conditioning of the matrix A , the problem may b e ill-p osed or ill-conditioned. In classical approaches, these issues are dealt with b y adjusting the original problem to a (sligh tly) mo dified optimization problem as discussed in Section 3.1 whose solution is meant to represent the original, as discussed in Section 3.2 . W e note that in the classical setting a solution to the inv erse problem is a v ector x as a result of solving an optimization problem. Suc h a solution is referred to as an p oint estimate b ecause it giv es one solution vector without providing any information ab out ho w reliable (or uncertain) the solution is [ 15 , 34 ]. On the other hand, the statistical in version approac h to inv erse problems pro vides an ensemble of solutions—defined b y the p osteriori distribution which not only p oint estimates can b e made but also their uncertain ty quantification [ 11 , 17 ]. Bayesian Optic al Flow with Unc ertainty Quantific ation 7 3.1. L e ast squar es solution The classical least squares solution to the inv erse problem is given by [ 12 ] x ` 2 = A † b , (20) where A † denotes the pseudo-in verse of A whic h can b e obtained from the singular v alue decomp osition of A [ 12 ]. Dep ending on the rank of A , the least squares solution x ` 2 is asso ciated with one of the minimization problems: ( min A x = b k x k 2 , if rank( A ) < n ; min x k A x − b k 2 , if rank( A ) = n . (21) Here k · k 2 denotes the ` 2 (Euclidean) norm. Let the true solution to Eq. ( 19 ) b e x ∗ , that is, b = A x ∗ + η . It follows that x ` 2 − x ∗ = A † η . (22) In practice, even when the matrix A has full column rank (rank( A ) = n ), the discrepancy b et ween the true and least squares solutions is typically dominated b y noise when some singular v alues of A are close to zero, rendering A an ill-conditioned matrix and the solution x ∗ unstable and sensitiv e to small changes in data [ 34 ]. 3.2. Tikhonov R e gularization A pow erful approach to resolv e the instabilities due to noise and the near-singularit y of A is to r e gularize the problem. In the classical Tikhono v regularization one adds a quadratic regularization term α x T L x for some prescrib ed matrix L to p enalize non- smo othness, giving rise to a regularized optimization problem [ 31 , 32 , 33 , 34 ]: min x  k A x − b k 2 2 + α x T L x  . (23) In the regularized problem the p ositive scalar parameter α con trols the weigh t of regularization and L is t ypically a symmetric p ositive definite matrix, b oth of which need to b e chosen appropriately for the problem to b e uniquely defined [ 31 , 32 , 33 , 34 ]. W e refer to the tw o terms k A x − b k 2 2 and α x > L x in ( 23 ) as data fidelity and solution r e gularity , resp ectiv ely . In a simplistic description, they can b e describ ed as “selecting” a solution x α that balances the desire to “solv e” A x = b and to b e “regular” as measured b y x > L x . The regularization parameter α therefore dictates the exten t to whic h the compromise is made b et w een the t wo. F or a given parameter α , we denote the corresp onding regularized solution b y x α = argmin x {k A x − b k 2 + α x > L x } . (24) By standard vector calculus, it can b e sho wn that x α is in fact a solution to the mo dified linear system ( A > A + αL ) x α = A > b , (25) Bayesian Optic al Flow with Unc ertainty Quantific ation 8 whic h is t ypically w ell-p osed for appropriate choices of L and α . When the matrices are large and sparse, Eq. ( 25 ) is often solv ed b y iterativ e metho ds rather than a direct matrix in version since the latter tends to b e numerically costly and unstable [ 12 ]. In Tikhonov regularization, a k ey issue is how to c ho ose the regularization parameter α appropriately . In theory , a “go o d” regularizer has the prop ert y that in the absence of noise, the solution to the regularized problem con verges to the true solution when the regularization parameter α → 0. How ev er, in practice, under the presence of noise, it is alwa ys a c hallenge to try to determine a go o d v alue for α . If α is to o small, the instability and sensitivit y of the original problem would still p ersist; whereas for to o large of an α the solution will b e o ver-regularized and not fit the data well. A go o d balance is th us paramount. Despite the existence and ongoing developmen t of man y comp eting metho ds for selecting α most of which fo cus on asymptotical optimalit y as the n um b er of data p oin ts approac h infinit y , none of them stands out as a best “natural” c hoice unless sp ecific priori information about the noise in the data are av ailable (see Chapter 7 of Ref. [ 34 ]). In the following section we discuss how the problem of selecting an exact regularization parameter is no longer required from a Ba yesian persp ective; instead, it suffices to start with some lo ose range of v alues represented by a probability distribution, unless more sp ecific knowledge is a v ailable ab out the problem in which case the distribution can b e chosen to reflect such information. 4. Statistical Inv ersion Approac h The statistical inv ersion approach to an inv erse problem starts by treating all v ariables as random v ariables (e.g., x , b and η in Eq. ( 19 )), and representing our knowledge (or absence of knowledge) of the unknowns as prior distributions. With observ ational data and a forward mo del (e.g., matrix A in Eq. ( 19 )), the “inv ersion” leads to an ensemble of p oints together with a p osterior distribution from whic h v arious statistical estimates can b e made [ 21 , 11 , 17 , 8 ]. The key in the inv ersion is to use the Bay es rule to express the p osterior distribution p ( x | b ), whic h is the conditional distribution of the “solution v ector” x giv en the observed data b , as [ 11 , 17 ] p ( x | b ) = 1 p ( b ) p ( b | x ) · p ( x ) . (26) Here the likeliho o d function p ( b | x ) is the probability density function (p df ) of the random v ariable b giv en x whic h is determined b y the underlying model; p ( x ) is the priori distribution of x ; and p ( b ) > 0 acts as a normalization constant which do es not affect the solution pro cedure or the final solution itself. Th us, in the statistical inv ersion formulation, each candidate solution x is asso ciated with the probability p ( x | b ) that is determined (up to a normalization constant 1 /p ( b )) once the lik eliho o d function and the prior distribution are giv en. F or a given in verse problem, the lik eliho o d function can be obtained using the underlying mo del suc h as Eq. ( 19 ) including the noise distribution. On the other hand, the prior distribution p ( x ) is typically constructed according to some prior kno wledge of the solution. Unlike Bayesian Optic al Flow with Unc ertainty Quantific ation 9 standard approaches which pro duce “p oint estimates”, in statistical inv ersion it is the p osterior distribution that is considered as the “solution” to the inv erse problem. Based on the p osterior distribution, one can further extract useful information such as p oint estimates and uncertain t y quan tification, by sampling from the distribution. Efficient sampling metho ds will b e review ed to w ard the end of this section. The unique feature of enabling information fusion and uncertaint y quantification has made the statistical inv ersion approac h to inv erse problems an attractive ven ue for the dev elopment of new theory and applications. In image pro cessing applications, it has b een utilized for many problems suc h as image denoising and deblurring [ 2 , 18 ], sparse signal reconstruction [ 23 ], and more recently attempted for optical flo w computation [ 9 ]. In particular, we note that our approach, although differen t in many of the tec hnical asp ects, shares a similar statistical inv ersion p ersp ectiv e as Ref. [ 9 ]. 4.1. Statistic al interpr etation of optic al flow obtaine d by Tikhonov r e gularization Under the statistical inv ersion framework, solution to an in verse problem is the p osterior distribution, from which p oint estimates can b e further obtained. Among these p oin t estimates, a particularly common one is the maximum a p osteriori (MAP) estimator, whic h is defined as x MAP = argmax x p ( x | p ) = argmin x {− ln p ( x | b ) } (27) = argmin x {− ln p ( b | x ) − ln p ( x ) } . (28) As noted in Refs. [ 2 , 17 ], the MAP estimator given b y Eq. ( 28 ) pro duces a v ector x that is iden tical to the solution of a Tikhonov regularization sp ecified in Eq. ( 24 ) up on appropriate choice of the mo del and prior p df. In particular, consider the mo del giv en b y Eq. ( 19 ) with indep endent and identically distributed (iid) Gaussian noise of v ariance λ − 1 . It follo ws that p ( b | x ) = p ( η ) ∝ exp  − λ 2 k A x − b k 2  , (29) where the sym b ol “ ∝ ” means “prop ortional to”. Under a Gaussian prior, p ( x ) ∝ exp  − δ 2 x > L x  , (30) the term − ln p ( x | b ) in the MAP estimator b ecomes − ln p ( x | b ) ∝ k A x − b k 2 + ( δ /λ ) x > L x . (31) The c hoice of α = δ /λ in the Tikhono v regularization then yields a solution x α that equals the vector x MAP giv en b y the same MAP estimator. Th us, with these assumptions of the form of the noise, the distribution of the prior, there is a logical bridge b et ween tw o different philosophies for in verse problems, in that a solution from Tikhono v regularization can b e interpreted as an MAP estimate of the p osteriori distribution under appropriate c hoices of prior distributions and likelihoo d functions. F urthermore, and crucially , as shown in the next section, more imp ortan t information exists in the statistical inv ersion framew ork. Sp ecifically , by sampling from Bayesian Optic al Flow with Unc ertainty Quantific ation 10 the p osterior distribution the statistical in v ersion approach allows for not only a p oint estimate but also other statistical prop erties asso ciated with the solution, in particular spread around a p oin t estimate which enables uncertaint y quan tification. 4.2. Choic e of Priors and Hyp erpriors In general, the statistical in version framew ork allo ws flexibilit y in c ho osing prior distributions ab out unkno wns and noise. When these prior distributions themselves con tain unknown parameters, these parameters can themselves b e though t of as random v ariables which follow hyperprior distributions. In this pap er we will consider some elemen tary c hoice of the priors and hyperpriors, with the main purp ose of showing how to setup the pro cedure for sampling the p osterior distribution. W e fo cus on (multiv ariate) Gaussian priors for b oth the unknown x and noise η . In the absence of knowledge of ho w “spread-out” the prior distributions are, w e use additional parameters for these priors, where these parameters are taken to b e random v ariables drawn from hyperprior distributions. This type of approach has b een previously adopted and implemen ted in image denoising applications [ 2 ]. In particular, w e consider a sp ecific class of the noise and prior distributions: ( lik eliho o d function: p ( b | x , λ ) ∝ λ m/ 2 exp  − λ 2 k A x − b k 2  , prior distribution: p ( x | δ ) ∝ δ n/ 2 exp  − δ 2 x > L x  . (32) Here the noise is assumed to b e additiv e, Gaussian, and indep enden t of the measured data, with v ariance λ − 1 , giving rise to the form of the likelihoo d function. On the other hand, the prior distribution is considered to b e Gaussian with co v ariance matrix ( δ L ) − 1 (matrix δ L is referred to as the pr e cision matrix). W e choose L = Q where Q is giv en b y Eq. ( 17 ) which corresp onds to a spatial regularization measure. As it turns out, this c hoice of L is closely related to the selection of prior according to a spatial Gaussian Mark ov random field which is common in tac kling spatial in verse problems [ 3 , 14 ]. As we p oin ted out, to completely sp ecify the p osterior distribution, we also need to c ho ose prior distributions for the parameters λ and δ . These are often called hyp erpriors . F ollowing Ref. [ 2 ], w e choose the priors for p ( λ ) and p ( δ ) to b e Gamma distributions, as ( p ( λ ) ∝ λ α λ − 1 exp( − β λ λ ) , p ( δ ) ∝ δ α δ − 1 exp( − β δ δ ) . (33) Suc h choice ensures that p ( λ ) and p ( δ ) are conjugate h yp er-priors, and mak es the design of sampling of p osterior distributions more conv enien t. In the absence of kno wledge of the v alues of λ and δ , one needs to choose the v alues of α λ , β λ , α δ , and β δ to ensure the distributions p ( λ ) and p ( δ ) to b e “wide”, allowing the Marko v chain to explore a p oten tially larger region of the parameter space. F ollowing Ref. [ 2 , 3 , 14 ], w e set α λ = α δ = 1 and β λ = β δ = 10 − 4 unless otherwise noted. W e tested other c hoice of parameters as w ell and they mainly affect the length of the transien t in the MCMC sampling pro cess and do not seem to ha ve a strong influence on the asymptotic outcome. Details of these will b e discussed in Section 5 along with n umerical examples. Bayesian Optic al Flow with Unc ertainty Quantific ation 11 4.3. Sampling fr om the Posterior Distribution In the statistical inv ersion formalism, once the form of the p osterior distribution is deriv ed, the remaining part of the w ork is devoted to efficient sampling from the p osterior distribution. Typically a Mark ov c hain Mon te Carlo (MCMC) sampling approach is adopted. The main idea is to generate a sequence of samples according to a prescrib ed Mark ov c hain whose unique stationary distribution is the desired p osterior distribution. Note that under the Gaussian priors and gamma hyperpriors as w e discussed ab o ve, the full conditional distributions that relate to the p osterior distribution are giv en by      p ( x | λ, δ, b ) ∝ exp  − λ 2 k A x − b k 2 − δ 2 x > L x  , p ( λ | x , δ, b ) ∝ λ m/ 2+ α λ − 1 exp  λ  − 1 2 k A x − b k 2 + β λ  , p ( δ | x , λ, b ) ∝ δ n/ 2+ α δ − 1 exp  δ  − 1 2 x > L x − β δ  . (34) In other w ords,      x | λ, δ, b ∼ N  ( λA > A + δ L ) − 1 λA > b , ( λA > A + δ L ) − 1  , λ | x , δ, b ∼ Γ  m/ 2 + α λ , 1 2 k A x − b k 2 + β λ  , δ | x , λ, b ∼ Γ  n/ 2 + α δ , 1 2 x > L x + β δ  . (35) W e adopt the blo ck Gibbs sampler developed in Refs. [ 2 , 17 ] as a sp ecific MCMC pro cedure to sample the p osterior distribution. In theory the sample distribution asymptotically conv erges to the true p osterior distribution. The approach con tains the follo wing steps. (i) Initialize δ 0 and λ 0 , and set k = 0. (ii) Sample x k ∼ N  ( λA > A + δ L ) − 1 λA > b , ( λA > A + δ L ) − 1  . (iii) Sample λ k +1 ∼ Γ  m/ 2 + α λ , 1 2 k A x k − b k 2 + β λ  . (iv) Sample δ k +1 ∼ Γ  n/ 2 + α δ , 1 2 ( x k ) > L ( x k ) + β δ  . (v) Set k ← k + 1 and return to Step (ii). Here the computational burden is mainly due to Step 1, whic h requires drawing samples from a m ultiv ariate Gaussian v ariable, which is equiv alent to solving the follo wing linear system at eac h iteration for x k : ( λ k A > A + δ k L ) x k = λ k A > b + w , where w ∼ N ( 0 , λ k A > A + δ k L ) . (36) F or large matrices, instead of a direct solv e using Gauss elimination, an iterativ e method is usually preferred. Among the v arious notable iterativ e methods suc h as Jacobi, Gauss- Seidel (G-S), and conjugate gradien t (CG) [ 22 ], w e adopted the CG for all the numerical exp erimen ts as rep orted in this pap er, with a starting vector of all zeros, maximum of 500 iterations, and error tolerance of 10 − 6 . Note that the computational b ottleneck in solving Eq. ( 36 ) can in principle b e tac kled using more efficien t methods, for example, b y exploring history of solutions with similar parameters to provide “go o d” initial guess or b y utilizing the Karhunen-Lo ` ev e expansion to reduce the dimensionality of the problem. Bayesian Optic al Flow with Unc ertainty Quantific ation 12 5. Results: Ba yesian Optical Flow from Statistical In version 5.1. Benchmark Flow Fields and Noisy Image Pairs T o b enc hmark the prop osed statistical inv ersion approach for Ba y esian optical flow, w e consider 5 qualitatively differen t flo w fields that span a broad div ersity of p ossibilities. The flow fields, which are all defined on the normalized spatial domain [ − 1 , 1] × [ − 1 , 1], are generally not divergence-free. In particular, the flow fields are defined as follo ws. Flo w field 1: h U ( x, y ) , V ( x, y ) i = h x, y i . Flo w field 2: h U ( x, y ) , V ( x, y ) i = h− y , x i . Flo w field 3: h U ( x, y ) , V ( x, y ) i = h y , sin( x ) i . Flo w field 4: h U ( x, y ) , V ( x, y ) i = h− π sin(0 . 5 π x ) cos(0 . 5 π y ) , π cos(0 . 5 π x ) sin(0 . 5 π y ) i . Flo w field 5: h U ( x, y ) , V ( x, y ) i = h− π sin( π x ) cos( π y ) , π cos( π x ) sin( π y ) i . W e first consider synthetic image pairs, where the first image is constructed b y the equation F ( x, y ) = 1 2 [cos( π x ) cos( π y ) + 1] . (37) Then, for eac h optical flow field, we generate the second image G using the equation g = f − f x u − f y v + η , (38) where f and g represent the v ectorization of F and G , resp ectiv ely , η denotes (m ultiv ariate) noise whose individual comp onen ts are indep enden tly drawn from a Gaussian distribution with zero mean and fixed standard deviation σ , and spatial deriv atives are n umerically implemen ted using the forward difference sc heme. All syn thetic images are constructed at the fixed resolution of 30-by-30 pixels, with uniform spacing in b oth directions, represented by a set of 30-by-30 matrices. 5.2. Bayesian Optic al Flow with Unc ertainty Quantific ation F or each image pair ( F , G ), w e adopt the MCMC-Gibbs sampling pro cedure and corresp onding c hoice of prior p df and hy erpriors presented in Section 4 to obtain an empirical p osterior distribution p ( U, V ) as the solution of the statistical inv ersion problem, which w e refer to as Bayesian optic al flow. Unlik e classical optical flow whic h pro vides a p oin t estimate, the Bay esian optical flo w can b e though t of as an ensemble of flow fields eac h asso ciated with some probability . F rom suc h an ensem ble and asso ciated probabilit y distribution (the p osterior distribution), w e can further extract useful information. F or example, the mean flo w field can b e computed using the sampling mean of the p osterior distribution, whic h is shown in Fig. 1 through Fig. 6 to compare with the true underlying flo w field. In particular, in each figure, the top ro w (a1-a3) sho ws the image data of the first image F (a1), and the second image G generated from Eq. ( 38 ) with no noise (a2) and with noise under standard deviation σ = 0 . 02 (a3), resp ectiv ely . The middle ro ws (b1-b3) show the true optical flow field (b1) compared with the inferred “mean” optical flow fields together with uncertaint y quan tification Bayesian Optic al Flow with Unc ertainty Quantific ation 13 from the MCMC samples (b2-b3). The uncertaint y regions are computed as follows. A t eac h p oin t z = ( x, y ), we construct a 2d normal p df N ( µ, Σ) b y using the sample mean µ and sample cov ariance Σ estimated from the MCMC samples after discarding the initial transien ts. This allo ws us to obtain a “mean” optical flow at p oin t z defined as h u ( z ) , v ( z ) i = h µ 1 , µ 2 i . Uncertaint y is quantified by computing a confidence region that con tains q probability mass of the fitted multiv ariate normal distribution given by [ 24 ] ( z − µ ) > Σ − 1 ( z − µ ) ≤ χ 2 2 ( q ) . (39) Here χ 2 2 ( q ) denotes the q -th quantile of the Chi-squared distribution with tw o degrees of freedom, that is, χ 2 2 ( q ) = K − 1 ( q ) where K is the cdf of χ 2 2 . These confidence regions (shaded ellipses) are shown for the zo omed-in plots for the inferred optical flow fields. Finally , the last row (c1-c3) in eac h figure sho ws ho w the MCMC pro cedure pro duces a distribution of the effective regularization parameter δ /λ . P anel (c1) sho ws the c hange of δ /λ o ver time in the MCMC sampling pro cedure, indicating conv ergence to a stationary distribution typically after a quic k initial transient. The remaining panels (c2) and (c3) sho w the distribution of δ /λ after discarding the initial transien t, for b oth the case of no noise (c2) and the case with noise (c3). W e p oin t out a few observ ations from the n umerical exp erimen ts. First, the estimated mean optical flo w compare reasonably well with the true flo w in all examples of the qualitatively different optical flo w fields, supporting the utility of the prop osed statistical in v ersion approac h [see panels (b1-b3) in all figures]. Secondly , we again point out that the MCMC pro cedure used in our statistical in version approach to optical flow do es not require an activ e prior c hoice of the regularization parameter. The MCMC samples seems to quic kly con verge to a stationary distribution for the effective parameter [panel (c1) in all figures], from which the distributions of parameters and solutions can b e determined. Finally , comparing to the noise free images, the estimation of optical flo w b ecomes less accurate when noise is added. It is worth mentioning that the statistical in version approac h in fact allows us to “predict” this difference without kno wing the ground-truth optical flo w, by quantifying and comparing the uncertain ty of solutions [panels (b2) and (b3) in all figures]. Note that although here w e only show the Ba yesian optical flow results based on Gaussian noise, we ha ve also p erformed sim ulations using uniform noise and Laplace-distributed noise, and found that the results are quite similar to what is shown in Fig. 1 to Fig. 6 . 5.3. Simulations on R e al Images Finally , w e conduct numerical exp erimen ts for the computation of Bay esian optical flo w from using real images together with the b enc hmark flow fields. In eac h example, w e consider matrix F defined by a real image from the Middlebury dataset ( http: //vision.middlebury.edu/stereo/ , also see [ 26 ]), resized to a fixed size of 60-by-60 pixels in grayscale, with in tensity normalized so that each pixel intensit y is in the range [0 , 1]. W e consider a total of six suc h images, as sho wn in Fig. 6 . Bayesian Optic al Flow with Unc ertainty Quantific ation 14 F or each real image F and a giv en flo w field h U, V i , w e generate a resulting second image G according to the flo w equation ( 38 ), where the noise is taken to b e a multiv ariate Gaussian distribution with zero mean and co v ariance matrix σ I with σ = 0 . 05. This c hoice of standard deviation ensures a non-negligible effect on the pixel intensities of the image pairs, since each F ij ∈ [0 , 1]. These “noisy second images” are sho wn in each of the first column of Fig. 7 to Fig. 12 . F or comparison, we also consider what we call a “true second image”, denoted as ¯ G , which is defined by the same flow equation ( 38 ) but in the absence of noise. These true second images are shown in the second column in Fig. 7 through Fig. 12 , for each one of the image example and flo w field. Thus, eac h baseline image F and flow field gives rise to a noisy image pair ( F , G ), and there is a total of 30 such pairs given the 6 real images and 5 flo w fields. F or each noisy image pair ( F , G ), w e use the same metho dology with the same c hoice of priors and h yp erpriors as in Sec. 5.2 to obtain a Ba y esian optical flow from sampling the posterior distribution–in practice because of the randomness of the MCMC pro cess, w e observ e that sometimes the sample distribution do es not “con verge” to a stationary distribution ev en after thousands of iterations; when this o ccurs we simply restart the MCMC from a differen t initial seed. Recall that different from a standard regularization approach which require careful choice of the regularization parameter, here suc h parameter is itself treated as a random v ariable that has its own prior which can b e taken to b e a “wide” distribution in the absence of additional knowledge. F rom the p osterior distribution, w e tak e the mean as an estimate of the a v erage flow field. T o test the usefulness of the reconstructed flow field, we use it together with the first image F to obtain an estimated second image ˆ G from equation ( 38 ), setting the noise term to b e zero. The estimated second image ˆ G is sho wn for eac h example in the third column of Fig. 7 through Fig. 12 . In terestingly , the estimated ˆ G seems to not only resemble the giv en noisy image G , but ev en more similar to the noiseless second image ¯ G . This observ ation is confirmed quan titatively , as sho wn in the last column of Fig. 7 through Fig. 12 , where the v alues of ˆ G ij are plotted against b oth those of G ij and of ¯ G ij . These results confirm the v alidit y of Ba yesian optical flow obtained by a statistical in v ersion approac h, and, additional suggests that accurate reconstruction of optical flow can b e p oten tially useful for image smo othing and denoising applications. 6. Discussion and Conclusions In this pap er we tak e a statistical inv ersion p ersp ectiv e to the optical flo w inference problem. F rom this p ersp ectiv e all relev an t v ariables in an otherwise standard in verse problem are treated as random v ariables, and the key is to form an efficient pro cess to construct and sample the p osterior distribution b y utilizing knowledge ab out the form of mo del, noise, and other prior information. F rom a Ba yesian p ersp ective, the v arious priors and the forward model combine to pro duce a p osterior distribution describing the propagation of prior information in con text of the problem. W e hav e shown that optical flow estimates given by traditional v ariational approaches suc h as the seminal Bayesian Optic al Flow with Unc ertainty Quantific ation 15 w ork developed by Horn and Sch unc k [ 16 ] can in fact b e in terpreted in the statistical in version framew ork under particular c hoices of mo dels, noise, and priors. Th us w e recap that there are ma jor adv an tages ov er the classical v ariational calculus approach to in verse problems where by necessity the ill-p osedness is dealt with by adding an ad ho c regularit y term that hop efully agrees with exp ected physical interpretation. F rom a Bay esian p ersp ectiv e, the ill-p osedness is dealt with naturally under the statistical in version framework by restating as a w ell-p osed extended problem in a larger space of probabilit y distributions [ 11 , 17 ]. This therefore naturally remo ves a k ey difficult y of ha ving to choose exactly an appropriate regularity parameter encoun tered in classical metho ds. Instead, in con trast to classical optical flow metho ds which only yield single solutions as “p oint estimates”, the statistical in version approach pro duces a distribution of p oints which can be sampled in terms of most appropriate estimators and also for uncertain ty quantification. Sp ecifically in the con text of an optic flow problem, w e exp ect a distribution of regularit y parameters, and corresp ondingly a distribution of optical flo w v ectors at eac h pixel. In this pap er we fo cused on a statistical inv ersion formulation with the c hoice of priors inspired from the classical Horn-Sch unc k functional. Although suc h c hoices app ear to b e reasonable for the synthetic flo w fields considered herein, real images and optical flow fields are muc h more complicated and, consequen tly , require additional efforts in forming the appropriate priors. Suc h priors can come from previous kno wledge of images and flow fields from similar systems, taken under similar scenes, or from other ph ysical measurements. In addition, differen t noise distributions should b e considered, and these will require mo difying the standard Horn-Sc hunc k framew ork, whic h assumed rigid b o dy motion and conserv ation of brightness. In particular, with regard to the regularization term, one example is to consider different p -norms other than the standard 2-norm, and p ossibly with a k ernel to w eigh in additional information ab out the physical em b edding of ob jects [ 27 ]. Other data fidelit y terms [ 5 , 19 , 6 ] can b e form ulated to describ e the ph ysics of the underlying application, for example for fluid and o ceanographic problems where a stream function, or ev en a quasi-static appro ximation to assume such ph ysics as corriolis can b e used; or in scenarios where div ergence-free flows are estimated by utilizing v orticit y-velocity formalism or more general data assimilation to ols [ 7 , 13 , 36 ]. Lik ewise, regularity in time and multiple time steps ma y b e appropriate [ 5 ], as these corresp ondingly more complex formulations nonetheless come bac k to a linear inv erse problem tenable in the framew ork of this pap er. Sp ecifically within the statistical in version framew ork dev elop ed in the current pap er, all of these could b e recast so that the data fidelit y term ma y be incorporated and should not require significan t mo dification of the general framework as in tro duced here, whic h w e plan in future work. Lik ewise, other numerical differen tiation and integration sc hemes can b e used as w ell in place of the simple forward difference used here. Bayesian Optic al Flow with Unc ertainty Quantific ation 16 References [1] G. Aub ert, R. Deriche and P . Kornprobst, Computing optical flow via v ariational techniques, SIAM J. Appl. Math. , 60 (1999), 156–182. [2] J. M. Bardsley , MCMC-based image reconstruction with uncertaint y quantification, SIAM J. Sci. Comput. , 34 (2012), A1316–A1332. [3] J. M. Bardsley , Gaussian Marko v random field priors for inv erse problems, Inverse Pr oblems and Imaging , 7 (2013), 397–416. [4] J. L. Barron, D. J. Fleet and S. S. Beauchemin, P erformance of optical flow tec hniques, International Journal of Computer Vision , 12 (1994), 43–77. [5] R. Basnay ake and E. M. Bollt, A Multi-Time Step Metho d to Compute Optical Flow with Scien tific Priors for Observ ations of a Fluidic System, in Er go dic The ory, Op en Dynamics, and Coher ent Structur es (eds. G. F ro yland and W. Bahsoun), Springer, 2014. [6] R. Basnay ake, A. Luttman and E. M. Bollt, A lagged diffusivity metho d for computing total v ariation regularized fluid flo w, Contemp. Math , 586 (2013), 57–64. [7] D. B´ er ´ eziat and I. Herlin, Solving ill-p osed image pro cessing problems using data assimilation, Numeric al A lgorithms , 56 (2011), 219–252. [8] L. Cav alier, Nonparametric statistical inv erse problems, Inverse Pr oblems , 24 (2008), 034004. [9] G. Chantas, T. Gk amas and C. Nikou, V ariational-Bay es optical flow, J. Math. Imaging Vision , 50 (2014), 199–213. [10] I. Cohen and I. Herlin, Optical flow and phase portrait metho ds for en vironmental satellite image sequences, in Pr o c. Eur op. Conf. Computer Vision , Cambridge, U.K., 1996, pp. 141–150. [11] A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. V ehtari and D. B. Rubin, Bayesian Data Anal ysis , 3 rd edition, Chapman & Hall/CRC, 2014. [12] G. H. Golub and C. F. V an Loan, Matrix Computations , 4 th edition, Johns Hopkins Universit y Press, 2012. [13] I. Herlin and D. B´ er´ eziat, Divergence-F ree Motion Estimation, Eur op e an Confer enc e on Computer Vision (ECCV 2012) , 15–27. [14] D. Higdon, A primer on space-time mo delling from a Ba yesian p erspective, in L os A lamos National L ab or atory, Statistic al Scienc es Gr oup, T e chnic al R ep ort , LA-UR-05-3097. [15] P . C. Hansen, R ank-Deficient and Discr ete Il l-Pose d Pr oblems: Numeric al Asp e cts of Line ar Inversion , SIAM, Philadelphia, 1998. [16] B. K. P . Horn and B. G. Sch unc k, Determining optical flow, Artificial Intel ligenc e , 17 (1981), 185–203. [17] J. Kaipio and E. Somersalo, Statistic al and Computational Inverse Pr oblems , Springer, 2005. [18] M. Lebrun, A. Buades and J. M. Morel, A nonlo cal Bay esian image denoising algorithm, SIAM J. Imaging Sci. , 6 (2013), 1665–1688. [19] A. Luttman, E. Bollt, R. Basnay ake and S. Kramer, A stream function approach to optical flow with applications to fluid transp ort dynamics, Pr o c. Appl. Math. Me ch. , 11 (2011), 855–856. [20] A. Luttman, E. M. Bollt, R. Basna yak e, S. Kramer and N. B. T ufillaro, A framework for estimating p oten tial fluid flow from digital imagery , Chaos , 23 (2013), 033134. [21] F. A. O’Sulliv an, A statistical p ersp ectiv e on ill-p osed inv erse problems, Statistic al Scienc e , 1 (1986), 502–527. [22] Y. Saad, Iter ative metho ds for sp arse line ar systems , 2 nd edition, SIAM, 2003. [23] M. W. Seeger and H. Nic kisc h, Large scale Bay esian inference and exp erimental design for sparse linear mo dels, SIAM J. Imaging Sci. , 1 (2011), 166–199. [24] M. Siotani, T olerance regions for a multiv ariate normal p opulation, Annals of the Institute of Statistic al Mathematics , 16 (1964), 135–153. [25] K. Souhila and A. Karim, Optical flo w based rob ot obstacle av oidance, International Journal of A dvanc e d R ob otic Systems , 4 (2007), 13–16. [26] D. Sun, S. Roth and M. J. Black, Secrets of optical flow estimation and their principles, Pr o c. Bayesian Optic al Flow with Unc ertainty Quantific ation 17 IEEE Conf. Computer Vision and Pattern R e c o gnition , (2010). [27] J. Sun, F. J Quevedo and E. Bollt, Data F usion Reconstruction of Spatially Embedded Complex Net w orks, arXiv pr eprint arXiv:1707.00731 , (2017). [28] R. Szeliski, Computer Vision , Springer 2011. [29] L. Le T arnec, F. Destremp es, G. Cloutier and D. Garcia, A pro of of conv ergence of the Horn- Sc h unc k optical flow algorithm in arbitrary dimension, SIAM J. Imaging Sci. , 7 (2014), 277–293. [30] A. M. Thompson, J. C. Brown, J. W. Kay and D. M. Titterington, A study of metho ds of c ho osing the smo othing parameter in image restoration by regularization, IEEE T r ans. Pattern Anal. Mach. Intel l. , 13 (1991), 326–339. [31] A. N. Tikhonov, Regularization of incorrectly p osed problems, Soviet Mathematics Doklady , 4 (1963), 1624–1627. [32] A. N. Tikhonov and V. Arsenin, Solutions of Il l-Pose d Pr oblems , Wiley , New Y ork, 1977. [33] A. N. Tikhonov, A. V. Goncharsky , V. V. Stepanov and A. G. Y agola, Numeric al Metho ds for the Solution of Il l-Pose d Pr oblems , Kluw er Academic Publishers, 1990. [34] C. R. V ogel, Computational Metho ds for Inverse Pr oblems , SIAM, Philadelphia, 2002. [35] J. W eick ert and C. Sc hn¨ orr, V ariational optic flo w computation with a spatio-temp oral smoothness constrain t, Journal of Mathematic al Imaging and Vision , 14 (2001), 245–255. [36] S. Zh uk, T. Tc hrakian, A. Akhriev, S. Lu and H. Hamann, Where computer vision can aid physics: dynamic cloud motion forecasting from satellite images, arXiv pr eprint arXiv:1710.00194 , (2017). Bayesian Optic al Flow with Unc ertainty Quantific ation 18 image 1 (F) (a1) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 true optical flow (b1) 0.1 0.3 0.5 0.7 0.1 0.3 0.5 0.7 image 2 (G) without noise (a2) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 inferred optical flow (b2) 0.1 0.3 0.5 0.7 0.1 0.3 0.5 0.7 1 2 3 δ / λ × 10 -4 0 200 400 600 800 1000 histogram of δ / λ (c2) 0 2000 4000 time step 10 -4 10 -2 10 0 10 2 time evolution of δ / λ (c1) (c1) without noise with noise image 2 (G) with noise (a3) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 inferred optical flow (b3) 0.1 0.3 0.5 0.7 0.1 0.3 0.5 0.7 0 2 4 6 8 δ / λ × 10 -3 0 200 400 600 800 1000 histogram of δ / λ (c3) Figure 1. Bay esian optical flow computed using statistical inv ersion, based on image pairs ( F , G ) under example flow field 1, where h U ( x, y ) , V ( x, y ) i = h x, y i . T op row (a1- a3): image data generated according to Eq. ( 38 ) for F (a1) and using Eq. ( 37 ) for G either without noise (a2) or with noise standard deviation σ = 0 . 02 (a3). Middle rows (b1-b3): true flow field (b1) and the mean flo w field from the p osterior distribution estimated using the MCMC samples. In the panels b elow (b2) and (b3), we also show the uncertain ty regions (as shaded ellipses), also from MCMC samples as giv en by Eq. ( 39 ). Bottom row (c1-c3): time evolution (c1) as w ell as the distribution (c2-c3) of the effective regularization parameter σ /λ , where the distributions are obtained after discarding the initial transient in the MCMC sampling pro cess. Bayesian Optic al Flow with Unc ertainty Quantific ation 19 image 1 (F) (a1) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 true optical flow (b1) 0.1 0.3 0.5 0.7 0.1 0.3 0.5 0.7 image 2 (G) without noise (a2) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 inferred optical flow (b2) 0.1 0.3 0.5 0.7 0.1 0.3 0.5 0.7 1 1.5 2 2.5 δ / λ × 10 -4 0 200 400 600 800 1000 histogram of δ / λ (c2) 0 2000 4000 time step 10 -6 10 -4 10 -2 10 0 10 2 time evolution of δ / λ (c1) (c1) without noise with noise image 2 (G) with noise (a3) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 inferred optical flow (b3) 0.1 0.3 0.5 0.7 0.1 0.3 0.5 0.7 0 2 4 6 8 δ / λ × 10 -3 0 200 400 600 800 1000 histogram of δ / λ (c3) Figure 2. Bay esian optical flow computed using statistical inv ersion, based on image pairs ( F , G ) under example flow field 2, where h U ( x, y ) , V ( x, y ) i = h− y , x i . T op row (a1-a3): image data generated according to Eq. ( 38 ) for F (a1) and using Eq. ( 37 ) for G either without noise (a2) or with noise standard deviation σ = 0 . 02 (a3). Middle rows (b1-b3): true flow field (b1) and the mean flo w field from the p osterior distribution estimated using the MCMC samples. In the panels b elow (b2) and (b3), we also show the uncertain ty regions (as shaded ellipses), also from MCMC samples as giv en by Eq. ( 39 ). Bottom row (c1-c3): time evolution (c1) as w ell as the distribution (c2-c3) of the effective regularization parameter σ /λ , where the distributions are obtained after discarding the initial transient in the MCMC sampling pro cess. Bayesian Optic al Flow with Unc ertainty Quantific ation 20 image 1 (F) (a1) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 true optical flow (b1) 0.1 0.3 0.5 0.7 0.1 0.3 0.5 0.7 image 2 (G) without noise (a2) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 inferred optical flow (b2) 0.1 0.3 0.5 0.7 0.1 0.3 0.5 0.7 1 2 3 δ / λ × 10 -4 0 200 400 600 800 1000 histogram of δ / λ (c2) 0 2000 4000 time step 10 -4 10 -2 10 0 10 2 time evolution of δ / λ (c1) (c1) without noise with noise image 2 (G) with noise (a3) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 inferred optical flow (b3) 0.1 0.3 0.5 0.7 0.1 0.3 0.5 0.7 0 2 4 6 8 δ / λ × 10 -3 0 200 400 600 800 1000 histogram of δ / λ (c3) Figure 3. Bay esian optical flow computed using statistical inv ersion, based on image pairs ( F , G ) under example flow field 3, where h U ( x, y ) , V ( x, y ) i = h y , sin( x ) i . T op row (a1-a3): image data generated according to Eq. ( 38 ) for F (a1) and using Eq. ( 37 ) for G either without noise (a2) or with noise standard deviation σ = 0 . 02 (a3). Middle rows (b1-b3): true flow field (b1) and the mean flo w field from the p osterior distribution estimated using the MCMC samples. In the panels b elow (b2) and (b3), we also show the uncertain ty regions (as shaded ellipses), also from MCMC samples as giv en by Eq. ( 39 ). Bottom row (c1-c3): time evolution (c1) as w ell as the distribution (c2-c3) of the effective regularization parameter σ /λ , where the distributions are obtained after discarding the initial transient in the MCMC sampling pro cess. Bayesian Optic al Flow with Unc ertainty Quantific ation 21 image 1 (F) (a1) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 true optical flow (b1) 0.1 0.3 0.5 0.7 0.1 0.3 0.5 0.7 image 2 (G) without noise (a2) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 inferred optical flow (b2) 0.1 0.3 0.5 0.7 0.1 0.3 0.5 0.7 2 4 6 8 δ / λ × 10 -5 0 200 400 600 800 1000 histogram of δ / λ (c2) 0 2000 4000 time step 10 -6 10 -4 10 -2 10 0 10 2 time evolution of δ / λ (c1) (c1) without noise with noise image 2 (G) with noise (a3) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 inferred optical flow (b3) 0.1 0.3 0.5 0.7 0.1 0.3 0.5 0.7 0 0.5 1 δ / λ × 10 -3 0 200 400 600 800 1000 histogram of δ / λ (c3) Figure 4. Ba y esian optical flow computed using statistical in version, based on image pairs ( F, G ) under example flow field 4, where h U ( x, y ) , V ( x, y ) i = h− π sin(0 . 5 π x ) cos(0 . 5 π y ) , π cos(0 . 5 π x ) sin(0 . 5 π y ) i . T op row (a1-a3): image data generated according to Eq. ( 38 ) for F (a1) and using Eq. ( 37 ) for G either without noise (a2) or with noise standard deviation σ = 0 . 02 (a3). Middle ro ws (b1-b3): true flo w field (b1) and the mean flow field from the p osterior distribution estimated using the MCMC samples. In the panels b elo w (b2) and (b3), we also show the uncertaint y regions (as shaded ellipses), also from MCMC samples as given by Eq. ( 39 ). Bottom ro w (c1-c3): time evolution (c1) as w ell as the distribution (c2-c3) of the effective regularization parameter σ /λ , where the distributions are obtained after discarding the initial transient in the MCMC sampling pro cess. Bayesian Optic al Flow with Unc ertainty Quantific ation 22 image 1 (F) (a1) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 true optical flow (b1) 0.1 0.3 0.5 0.7 0.1 0.3 0.5 0.7 image 2 (G) without noise (a2) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 inferred optical flow (b2) 0.1 0.3 0.5 0.7 0.1 0.3 0.5 0.7 0 2 4 6 δ / λ × 10 -5 0 2000 4000 6000 8000 10000 histogram of δ / λ (c2) 0 1 2 3 time step × 10 4 10 -5 10 0 10 5 time evolution of δ / λ (c1) (c1) without noise with noise image 2 (G) with noise (a3) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 inferred optical flow (b3) 0.1 0.3 0.5 0.7 0.1 0.3 0.5 0.7 0 2 4 δ / λ × 10 -4 0 2000 4000 6000 8000 10000 histogram of δ / λ (c3) Figure 5. Ba y esian optical flow computed using statistical in version, based on image pairs ( F, G ) under example flow field 5, where h U ( x, y ) , V ( x, y ) i = h− π sin( π x ) cos( π y ) , π cos( π x ) sin( π y ) i . T op ro w (a1-a3): image data generated according to Eq. ( 38 ) for F (a1) and using Eq. ( 37 ) for G either without noise (a2) or with noise standard deviation σ = 0 . 02 (a3). Middle rows (b1-b3): true flow field (b1) and the mean flow field from the p osterior distribution estimated using the MCMC samples. In the panels b elow (b2) and (b3), w e also sho w the uncertain t y regions (as shaded ellipses), also from MCMC samples as given by Eq. ( 39 ). Bottom ro w (c1-c3): time evolution (c1) as w ell as the distribution (c2-c3) of the effective regularization parameter σ /λ , where the distributions are obtained after discarding the initial transient in the MCMC sampling pro cess. Bayesian Optic al Flow with Unc ertainty Quantific ation 23 (a) (b) (c) (d) (e) (f) Figure 6. Real images used for the computations of Bay esian optical flo w in Sec. 5.3 . Bayesian Optic al Flow with Unc ertainty Quantific ation 24 noisy second image, G flow field 1 true second image, ¯ G estimated second image, ˆ G 0 0.5 1 pixe l inten sities, ˆ G ij 0 0.2 0.4 0.6 0.8 1 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 2 true second image, ¯ G estimated second image, ˆ G 0 0.5 1 pixe l inten sities, ˆ G ij 0 0.2 0.4 0.6 0.8 1 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 3 true second image, ¯ G estimated second image, ˆ G 0 0.5 1 pixe l inten sities, ˆ G ij 0 0.2 0.4 0.6 0.8 1 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 4 true second image, ¯ G estimated second image, ˆ G -0.5 0 0.5 1 1.5 pixe l inten sities, ˆ G ij -0.5 0 0.5 1 1.5 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 5 true second image, ¯ G estimated second image, ˆ G -0.5 0 0.5 1 pixe l inten sities, ˆ G ij -0.5 0 0.5 1 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij Figure 7. Optical flow computed from noisy image pairs ( F , G ), where F is given by the image shown in Fig. 6 (a). Each row corresp onds to the choice of a differen t flow field, each of which is defined in Sec. 5.1 . F rom left to right in each ro w: noisy second image G obtained from the flow equation ( 38 ) using Gaussian noise with σ = 0 . 05, true second image ¯ G from the same flow equation in the absence of noise, estimated ˆ G based on using the mean flow field, and plots of ˆ G against G versus ˆ G against ¯ G . Bayesian Optic al Flow with Unc ertainty Quantific ation 25 noisy second image, G flow field 1 true second image, ¯ G estimated second image, ˆ G 0 0.5 1 pixe l inten sities, ˆ G ij 0 0.2 0.4 0.6 0.8 1 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 2 true second image, ¯ G estimated second image, ˆ G 0 0.5 1 pixe l inten sities, ˆ G ij 0 0.5 1 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 3 true second image, ¯ G estimated second image, ˆ G 0 0.5 1 pixe l inten sities, ˆ G ij 0 0.2 0.4 0.6 0.8 1 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 4 true second image, ¯ G estimated second image, ˆ G -1 0 1 pixe l inten sities, ˆ G ij -1 -0.5 0 0.5 1 1.5 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 5 true second image, ¯ G estimated second image, ˆ G -1 0 1 pixe l inten sities, ˆ G ij -1 -0.5 0 0.5 1 1.5 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij Figure 8. Optical flow computed from noisy image pairs ( F , G ), where F is given by the image shown in Fig. 6 (b). Each row corresp onds to the choice of a differen t flow field, each of which is defined in Sec. 5.1 . F rom left to right in each ro w: noisy second image G obtained from the flow equation ( 38 ) using Gaussian noise with σ = 0 . 05, true second image ¯ G from the same flow equation in the absence of noise, estimated ˆ G based on using the mean flow field, and plots of ˆ G against G versus ˆ G against ¯ G . Bayesian Optic al Flow with Unc ertainty Quantific ation 26 noisy second image, G flow field 1 true second image, ¯ G estimated second i mage, ˆ G 0 0.5 1 pixe l inten sities, ˆ G ij 0 0.5 1 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 2 true second image, ¯ G estimated second i mage, ˆ G 0 0.5 1 pixe l inten sities, ˆ G ij 0 0.5 1 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 3 true second image, ¯ G estimated second i mage, ˆ G 0 0.5 1 pixe l inten sities, ˆ G ij 0 0.5 1 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 4 true second image, ¯ G estimated second i mage, ˆ G -0.5 0 0.5 1 1.5 pixe l inten sities, ˆ G ij -0.5 0 0.5 1 1.5 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 5 true second image, ¯ G estimated second i mage, ˆ G -0.5 0 0.5 1 1.5 pixe l inten sities, ˆ G ij -0.5 0 0.5 1 1.5 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij Figure 9. Optical flow computed from noisy image pairs ( F , G ), where F is given by the image shown in Fig. 6 (c). Eac h row corresp onds to the choice of a differen t flow field, each of which is defined in Sec. 5.1 . F rom left to right in each ro w: noisy second image G obtained from the flow equation ( 38 ) using Gaussian noise with σ = 0 . 05, true second image ¯ G from the same flow equation in the absence of noise, estimated ˆ G based on using the mean flow field, and plots of ˆ G against G versus ˆ G against ¯ G . Bayesian Optic al Flow with Unc ertainty Quantific ation 27 noisy second image, G flow field 1 true second image, ¯ G estimated second image, ˆ G 0 0.5 1 pixe l inten sities, ˆ G ij 0 0.2 0.4 0.6 0.8 1 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 2 true second image, ¯ G estimated second image, ˆ G 0 0.5 1 pixe l inten sities, ˆ G ij 0 0.5 1 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 3 true second image, ¯ G estimated second image, ˆ G 0 0.5 1 pixe l inten sities, ˆ G ij 0 0.5 1 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 4 true second image, ¯ G estimated second image, ˆ G 0 1 2 pixe l inten sities, ˆ G ij -0.5 0 0.5 1 1.5 2 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 5 true second image, ¯ G estimated second image, ˆ G -1 0 1 2 pixe l inten sities, ˆ G ij -1 0 1 2 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij Figure 10. Optical flow computed from noisy image pairs ( F , G ), where F is given b y the image shown in Fig. 6 (d). Each row corresp onds to the choice of a differen t flow field, each of which is defined in Sec. 5.1 . F rom left to right in each ro w: noisy second image G obtained from the flow equation ( 38 ) using Gaussian noise with σ = 0 . 05, true second image ¯ G from the same flow equation in the absence of noise, estimated ˆ G based on using the mean flow field, and plots of ˆ G against G versus ˆ G against ¯ G . Bayesian Optic al Flow with Unc ertainty Quantific ation 28 noisy second image, G flow field 1 true second image, ¯ G estimated second i mage, ˆ G 0 0.5 1 pixe l inten sities, ˆ G ij -0.2 0 0.2 0.4 0.6 0.8 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 2 true second image, ¯ G estimated second i mage, ˆ G 0 0.5 1 pixe l inten sities, ˆ G ij -0.2 0 0.2 0.4 0.6 0.8 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 3 true second image, ¯ G estimated second i mage, ˆ G 0 0.5 1 pixe l inten sities, ˆ G ij -0.2 0 0.2 0.4 0.6 0.8 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 4 true second image, ¯ G estimated second i mage, ˆ G -0.5 0 0.5 1 pixe l inten sities, ˆ G ij -0.5 0 0.5 1 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 5 true second image, ¯ G estimated second i mage, ˆ G -0.5 0 0.5 1 pixe l inten sities, ˆ G ij -0.5 0 0.5 1 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij Figure 11. Optical flow computed from noisy image pairs ( F , G ), where F is given b y the image shown in Fig. 6 (e). Eac h row corresp onds to the choice of a differen t flow field, each of which is defined in Sec. 5.1 . F rom left to right in each ro w: noisy second image G obtained from the flow equation ( 38 ) using Gaussian noise with σ = 0 . 05, true second image ¯ G from the same flow equation in the absence of noise, estimated ˆ G based on using the mean flow field, and plots of ˆ G against G versus ˆ G against ¯ G . Bayesian Optic al Flow with Unc ertainty Quantific ation 29 noisy second image, G flow field 1 true second image, ¯ G estimated second image, ˆ G 0 0.5 1 pixe l inten sities, ˆ G ij -0.2 0 0.2 0.4 0.6 0.8 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 2 true second image, ¯ G estimated second image, ˆ G 0 0.5 1 1.5 pixe l inten sities, ˆ G ij 0 0.5 1 1.5 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 3 true second image, ¯ G estimated second image, ˆ G 0 0.5 1 pixe l inten sities, ˆ G ij 0 0.2 0.4 0.6 0.8 1 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 4 true second image, ¯ G estimated second image, ˆ G -1 0 1 pixe l inten sities, ˆ G ij -1 -0.5 0 0.5 1 1.5 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij noisy second image, G flow field 5 true second image, ¯ G estimated second image, ˆ G -0.5 0 0.5 1 1.5 pixe l inten sities, ˆ G ij -0.5 0 0.5 1 1.5 ˆ G ij vs. ¯ G ij ˆ G ij vs. G ij Figure 12. Optical flow computed from noisy image pairs ( F , G ), where F is given b y the image shown in Fig. 6 (f ). Eac h row corresponds to the choice of a different flo w field, each of which is defined in Sec. 5.1 . F rom left to right in each ro w: noisy second image G obtained from the flow equation ( 38 ) using Gaussian noise with σ = 0 . 05, true second image ¯ G from the same flow equation in the absence of noise, estimated ˆ G based on using the mean flow field, and plots of ˆ G against G versus ˆ G against ¯ G .