# Scene Flow Estimation Based on Adaptive Anisotropic Total Variation Flow-Driven Method

Authors: Xuezhi Xiang; Rongfang Zhang; Mingliang Zhai; Deguang Xiao; Erwei Bai
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Scene flow estimation based on disparity and optical flow is a challenging task. We present a novel method based on adaptive anisotropic total variation flow-driven method for scene flow estimation from a calibrated stereo image sequence. The basic idea is that diffusion of flow field in different directions has different rates, which can be used to calculate total variation and anisotropic diffusion automatically. Brightness consistency and gradient consistency constraint are employed to establish the data term, and adaptive anisotropic flow-driven penalty constraint is employed to establish the smoothness term. Similar to the optical flow estimation, there are also large displacement problems in the estimation of the scene flow, which is solved by introducing a hierarchical computing optimization. The proposed method is verified by using the synthetic dataset and the real scene image sequences. The experimental results show the effectiveness of the proposed algorithm.


---

## Body

## 1. Introduction

Scene flow was first introduced by Vedula et al. [1] in 1999. It is defined as a dense 3D motion field that describes the change of structured light in the image surfaces, which can be estimated from a calibrated stereo image sequence. We can obtain stereo image sequence by using binocular stereo vision sensor, like BumblebeeX3. As we all know, optical flow can be seen as the projection of the moving object on the retina or the camera’s sensor. So scene flow can be regarded as a 3D extension of optical flow. Reliable 3D motion flow field can be used in many applications. For example, Menze and Geiger [2] focus on the autonomous driving with 3D scene flow estimation. Lenz et al. [3] used scene flow for moving object detection in urban environments. Waizenegger et al. [4] present a real-time approach for 3D structure estimation of human bodies. In general, scene flow estimation techniques can be divided into two categories. One is variational minimization method of scene flow based on the optical flow and the disparity from stereo image sequence [5]. The other is 3D point cloud parametrization method based on the 3D structure to estimate the scene flow, which allows us to directly estimate the desired unknowns [6]. Under these two frameworks, a lot of researches on optical flow estimation can be used for scene flow estimation.In the estimation of optical flow, many regularization approaches have been derived to solve the ill-posed problem. These methods can be done by minimizing an energy function(1)min∫ΩΨdataIx·u+Iy·v+It+Ψsmoothu,vdxdy,where Ψsmooth(u,v) is a regularization term which penalizes variations u and v. Ω refers to the domain of the image plane. Ψdata and Ψsmooth are convex regularization functions and used to guarantee convexity and differentiability of the energy function.Many smoothing regularizations are characterized as below. Homogenous regularization is named by Bruhn et al. [7], because the approach is a linear diffusion equation which is defined by Horn and Schunck [8]. Inspired by image filtering, Alvarez and Esclarin [9] use the squared magnitude of the image gradient to control the smoothness. Perona and Malik [10] propose an anisotropic smoothness method, which can preserve flow edges for a large image gradient and reduce the smooth strength in the consistency area. Nagel and Enkelmann [11] only smooth flow along the image edge and reduce the smooth strength in the direction perpendicular to image edge. In other words, this is the “oriented smoothness” method, which is driven by image gradient. Weickert and Schnörr [12] propose a systematic classification of rotation invariant convex regularizers, which provides a unifying framework for image-driven and flow-driven methods, as well as isotropic and anisotropic regularizers.Inspired by optical flow regularizers, we propose a scene flow estimation method, which can take the advantages of anisotropic total variation and adaptive regularizer. Different from the structure tensor, our anisotropic diffusion applies robust penalty function to the component of gradient terms in different direction. A weight function used by Blomgren and Chan [13] for image segmentation is introduced to reduce diffusion of flow field in the large gradient area.


## 2. Related Work

Binocular camera can be used to build the stereo vision system that is used to obtain the stereo image sequences. Li and Sclaroff [14] first estimated scene flow under binocular setting as a joint framework including optical flow estimation and stereo matching. Huguet and Devernay [5] proposed a common variational framework for scene flow estimation by coupling optical flow and stereo matching. A rotating sphere dataset with ground truth is provided publicly in this paper, which allows us to quantitatively evaluate the accuracy of the estimation. After that, KITTI benchmark was published for autonomous driving [2], and scene flow estimation can also be evaluated using this dataset. 400 dynamic scenes were annotated with detailed 3D CAD models and revealed great challenges for scene flow algorithms. Wedel et al. [15] modified Huguet and Devernay’s method [5] by decoupling the stereo and motion estimation. Optical flow consistency in both views and the stereo consistency in time t+1 are taken as data term, while the smoothness term of disparity change was separated from the regularization term. With GPU’s development, they achieved the promising accuracy results with frame-rate at 5 fps. Basha et al. [6] proposed a 3D point cloud parametrization method which enforces multiview geometric consistency based on brightness consistency and piecewise smoothness assumption. Menze and Geiger [2] used a novel model and dataset for scene flow estimation, which represents each element in the scene by its rigid motion parameters. And each element can be achieved by super pixel segmentation. Vogel et al. [16] proposed a piecewise rigid model, which combines rigidly moving 3D planes constraint and segmentation regularization, and their method yields better results.


## 3. Proposed Method of Scene Flow Estimation

### 3.1. Variational Scene Flow Framework

Our goal is unify the data term and smooth term in a variational framework and obtain the scene flow and 3D structure. Our method follows the variational formulation proposed by Huguet and Devernay [5]. It has the following common form:(2)Eu,v,d0,dt=∫ΩEdata+αEsmoothdxdy.


### 3.2. Data Term

We first obtain a set of stereo image sequences, whose disparity is only in the horizontal direction. Two image sequences between two time steps are shown in Figure1. According to the time and space constraint from 4 images, data term can be divided into 4 parts(3)Edata=∫ΩOlEl+OrEr+Od0Ed0+OdtEdtdxdy,where Ω usually refers to the domain of the image plane. O is an occlusion factor, which equals 0 when the pixel is occluded and equals 1 otherwise. Each part contains brightness constancy assumption and gradient constancy assumption to overcome the nonlinear illumination problem.Figure 1
Data term in the two-view case. Motion of a single pixel point in two consecutive frames of the image sequences. Brightness and gradient can be considered equal approximately in every image.In other words, the brightness and the gradient of the pixel are nearly constant in the short time and in different views. According to the corresponding relationship, brightness and gradient constancy assumption are given by(4)Ilx=Ilx+w=Irx+d0+Irx+w+dt,∇Ilx=∇Ilx+w=∇Irx+d0+∇Irx+w+dt,where ∇ is the vector of gradient. We introduce the following abbreviation to explain these 4 terms.(5)ΔI,x;I′,x+dx=I′x+dx-Ix2+γ∇I′x+dx-∇Ix2.The 4 components of the data term can be written as(6)El=ΨΔIl,x;Il,x+w,Er=ΨΔIr,x+d0;Ir,x+w+dt,Ed0=ΨΔIl,x;Ir,x+d0,Edt=ΨΔIl,x+w;Ir,x+w+dt,where El and Er are the left and the right data term energy function of optical flow, respectively. Ed0 is the disparity term that represents the stereo matching at time t0. Edt is the disparity term at time t+1. Penalty function (Ψ(s2)=s2+ε2,ε=0.001) is introduced to guarantee convexity and differentiability of the whole data term energy function and remove the outliers at the same time. Outliers are the pixels that are not consistent with the pixels affected by noise, illumination change, and occlusion. γ can be determined according to the intensity of illumination change or just rely on the experiences.The most obvious thing should be noted here is that if there is an error in any term ofu,v,d0, and dt, the error would be propagated to others. We minimize the 4 data terms together to prevent the error from this situation.


### 3.3. Smooth Term Based on Anisotropic Total Variation

Huguet and Devernay [5] used flow-driven regularization that uses a robust and convex function to keep convexity and differentiability of the smooth term(7)Esmoothu,v,d0,dt=∫ΩΨ∇u2+∇v2+ρ∇d0-dt2+μ∇d02dxdy.Here Ψ is used to remove the outliers and preserve the discontinuities of the flow field. The magnitude of the gradient is(8)Ψ∇u=ux2+uy2+ε2.We further derive its scene flow form by minimizing the following problems. In order to simplify the equation, we only consider the 2D case(9)F=∫ΩΨ∇u2+∇v2dxdy=∫Ωux2+uy2+vx2+vy2+ε2dxdy.It corresponds to the diffusion equation as(10)∂F∂u=divΨ′∇u2+∇v2∇u=∂∂xuxux2+uy2+vx2+vy2+ε2+∂∂yuyux2+uy2+vx2+vy2+ε2,∂F∂v=divΨ′∇u2+∇v2∇v=∂∂xvxux2+uy2+vx2+vy2+ε2+∂∂yvyux2+uy2+vx2+vy2+ε2.In the directions of u and v, the weights are the same and are related to the sum of the squares of the gradient. We propose a new scheme that regularization includes information about the direction of the flow and we define it as anisotropic total variation smoothness term. In the 2D case,(11)TVanisotropic=∑i,j∇u+∑i,j∇v=∑i=1M∑j=1Nui,j-ui,j+1+ui,j-ui+1,j+∑i=1M∑j=1Nvi,j-vi,j+1+vi,j-vi+1,j.But for outliers, robustness function still should be used. Here we show that it is anisotropic. Similarly, in the 2D case, we define the following equation:(12)F=∫ΩΨux2+Ψuy2+Ψvx2+Ψvy2dxdy=∫Ωux2+ε2+uy2+ε2+vx2+ε2+vy2+ε2dxdy.We minimize it and obtain the following equation:(13)∂F∂u=∂∂xuxux2+ε2+∂∂yuyuy2+ε2,∂F∂v=∂∂xvxvx2+ε2+∂∂yvyvy2+ε2.Obviously, the weights only depend on the direction of each gradient of ux or vx and are not related to other direction, which can be regarded as anisotropic diffusion of flow field. In this way, we can control the rate of diffusion according to the gradient of flow field so that we can preserve the sharp edge at the position of nonflat area. And in the flat region, it has the same strength of smooth term. Then a new smooth term can be defined as(14)Esmoothu,v,d0,dt=∫ΩΨux2+Ψuy2+Ψvx2+Ψvy2+ρΨd0-dtx2+Ψd0-dty2+μΨd0x2+Ψd0y2dxdy,where ρ and μ are empirical parameters.


### 3.4. Smooth Term Based on Adaptive Anisotropic Total Variation Flow-Driven Method

A major deficiency of the total variation that Rudin et al. [17] pointed out is staircase in the steady state solution. Marquina and Osher [18] proposed a modified method in 2000, but they just relieved this phenomenon. Different from total variation, there is no such problem with L2 norm. Blomgren and Chan [13] used the following norm in the image restoration:(15)Rpu=∫Ω∇up∇udxdy,where p function monotonously decreases and meets the following requirement:(16)p∇u⟶2,∇u⟶01,∇u⟶∞.At the edge position, we have(17)Rpu≈∫Ω∇udxdy.In the flat region, we have(18)Rpu≈∫Ω∇u2dxdy.In other words, this norm can choose from L1 to L2 according to the gradient magnitude automatically. To achieve this, a rational form can be chosen as follows:(19)p∇u=∇u+2∇u+1.The Euler-Lagrange equation for (15) is(20)∂R∂u=divp∇u·∇up∇u-1+ln⁡p∇up′∇u·∇up∇u∇u∇u.Define a function to represent the equation of the bracket in (20):(21)q∇u=p∇u·∇up∇u-1+ln⁡p∇up′∇u·∇up∇u.Then (20) can be simplified as(22)∂R∂u=divq∇u∇u∇u.Equation (23) can be derived from (21):(23)q∇u=2∇u,∇u⟶01,∇u⟶∞.In the smooth region of flow, the regularization is similar to linear smoothing, and total variation model will be used closing to the flow edge. Combining the robust function Ψ(s2)=s2+ε2, our scene flow smooth term can be written as(24)Esmoothu,v,d0,dt=∫ΩΨux2pux+Ψuy2puy+Ψvx2pvx+Ψvy2pvy+ρΨd0-dtx2pd0-dtx+Ψd0-dty2pd0-dty+μΨd0x2pd0x+Ψd0y2pd0ydxdy.


### 3.5. Optimization

#### 3.5.1. Euler-Lagrange Equations

Before the introduction of minimizing the energy function, we first introduce the following abbreviations for the left image. The similar abbreviation is also applicable to the right image.(25)Ilx≔∂xIlx+w,Ilxz≔∂xIlx+w-∂xIlx,Ily≔∂yIlx+w,Ilyz≔∂yIlx+w-∂yIlx,Ilz≔Ilx+w-Ilx,Ilyy≔∂yy2Ilx+w,Ilxx≔∂xx2Ilx+w,Ilxy≔∂xy2Ilx+w,Ilt+1≔Ilx+w.We also use abbreviations to replace the partial derivatives of the data term.(26)Ψ′fl≔∂xΨIlx+w-Ilx2+γ∇Ilx+w-∇Ilx2,Ψ′fr≔∂xΨIrx+w+dt-Irx+d02+γ∇Irx+w+dt-∇Irx+d02,Ψ′st≔∂xΨIrx+w+dt-Ilx+w2+γ∇Irx+w+dt-∇Ilx+w2.It should be noted that the regularization term we used is anisotropic, so the Euler-Lagrange equations are not the same as (21). Considering anisotropic total variation of individual variables(27)F=Ψux2pux+Ψuy2puy,we can get the following Euler equations:(28)∂F∂ux=2ux·pux·ψ′ux2·ψux2pux-1+ln⁡ψux2·p′ux·Ψux2puxuxux,∂F∂uy=2uy·puy·Ψ′uy2·Ψuy2puy-1+ln⁡Ψuy2·p′uy·Ψuy2puyuyuy.We define the q1 function as follows:(29)q1s=2s·ps·Ψ′s2·Ψs2ps-1+ln⁡Ψs2·p′s·Ψs2ps,and then (28) can be simplified as(30)∂F∂ux=q1uxuxux,∂F∂uy=q1uyuyuy.According to the variational principle, we minimize the energy function and achieve(31)∇Eu,v,d0,dt=0.Equation (31) corresponds to the 4 Euler equations. We calculate the partial derivation of E with respect to u and let it be equal to 0 to get the following equation:(32)βlΨ′fl·IlxIlz+γIlxxIlxz+IlxyIlyz+βrΨ′fr·IrxIrz+γIrxxIrxz+IrxyIryz+βstΨ′st·Irt+1-Ilt+1Irx-Ilx+γIrx-IlxIrxx-Ilxx+Iry-IlyIrxy-Ilxy-α∂∂xq1uxuxux+∂∂yq1uyuyuy=0.The partial derivatives of E with respect to v,d0-dt, and d0 produce similar equations as given below:(33)βlΨ′fl·IlyIlz+γIlxyIlxz+IlyyIlyz+βrΨ′fr·IryIrz+γIrxyIrxz+IryyIryz+βstΨ′st·Irt+1-Ilt+1Iry-Ily+γIrx-IlxIrxy-Ilxy+Iry-IlyIryy-Ilyy-α∂∂xq1vxvxvx+∂∂yq1vyvyvy=0,βrΨ′fr·IrxIrz+γIrxxIrxz+IrxyIryz+βstΨ′st·IrxIrt-Ilt+γIrxxIrx-Ilx+IrxyIry-Ily-αρ∂∂xq1dt-d0xdt-d0xdt-d0x+αρ∂∂yq1dt-d0ydt-d0ydt-d0y=0,βrΨ′fr·IrxtIrz+γIrxxtIrxz+IrxytIryz+βstΨ′st·IrxtIrz+γIrxxtIrxt-Ilxt+IrxytIryt-Ilyt-αρ+μ∂∂xq1d0xd0xd0x+∂∂yq1d0yd0yd0y=0.


#### 3.5.2. Multiresolution Strategy

In order to solve the problem of large displacements in the scene and avoid local minima, we introduce the multiresolution method. Referring to Brox et al. [19], the sampling factor is set to 0.9 to avoid a larger span of adjacent layers and local minima. At the same time, we also have a serious problem in the scene flow estimation. If we start to calculate from low resolution layer, there will be a false disparity which is passed to the high resolution layer. For this reason, we start to calculate scene flow from middle resolution layer. A schematic diagram that uses the pyramid hierarchy is given in Figure 2.Figure 2
Multiresolution strategy for scene flow computation. The left side pyramid is the image at timet. The right side pyramid is the image at time t+1.In each level of the pyramid, we use the warping method. For the two frames of theK level, I0k is obtained according to a certain sampling standard. Itk can be obtained by motion compensation based on the image Itk-1 in upper level. We only need to minimize the residual of each layer of the scene flow and the error generated by the large displacement can be avoided. The evaluation of scene flow requires minimizing partial differential equations (32)-(33). This is done by using Euler-Lagrange equations. We use a nested fixed point iteration scheme suggested by Brox et al. [19] and linearize equations in each layer of pyramid. The inner loop is used to linearize brightness constancy assumption and outer loop linearizes the robust penalty function Ψ. We use SOR (Super Relaxation Iteration) as an iterative method. The full algorithm for scene flow estimation is given in the algorithm scheme (see Algorithm 1).Algorithm 1:Scene flow estimation.
Ensure: The images at time t and t+1 are in the same level Require: a>b≥c≥0 for level = c to b  do Left flow (ul,vl) adaptive anisotropic flow-driven Right flow (ur,vr) adaptive anisotropic flow-driven end for stereo (d0,dt) for level = b to a  do scene flow (ul,vl,d0,dt) adaptive anisotropic total variation flow-driven end for


## 3.1. Variational Scene Flow Framework

Our goal is unify the data term and smooth term in a variational framework and obtain the scene flow and 3D structure. Our method follows the variational formulation proposed by Huguet and Devernay [5]. It has the following common form:(2)Eu,v,d0,dt=∫ΩEdata+αEsmoothdxdy.


## 3.2. Data Term

We first obtain a set of stereo image sequences, whose disparity is only in the horizontal direction. Two image sequences between two time steps are shown in Figure1. According to the time and space constraint from 4 images, data term can be divided into 4 parts(3)Edata=∫ΩOlEl+OrEr+Od0Ed0+OdtEdtdxdy,where Ω usually refers to the domain of the image plane. O is an occlusion factor, which equals 0 when the pixel is occluded and equals 1 otherwise. Each part contains brightness constancy assumption and gradient constancy assumption to overcome the nonlinear illumination problem.Figure 1
Data term in the two-view case. Motion of a single pixel point in two consecutive frames of the image sequences. Brightness and gradient can be considered equal approximately in every image.In other words, the brightness and the gradient of the pixel are nearly constant in the short time and in different views. According to the corresponding relationship, brightness and gradient constancy assumption are given by(4)Ilx=Ilx+w=Irx+d0+Irx+w+dt,∇Ilx=∇Ilx+w=∇Irx+d0+∇Irx+w+dt,where ∇ is the vector of gradient. We introduce the following abbreviation to explain these 4 terms.(5)ΔI,x;I′,x+dx=I′x+dx-Ix2+γ∇I′x+dx-∇Ix2.The 4 components of the data term can be written as(6)El=ΨΔIl,x;Il,x+w,Er=ΨΔIr,x+d0;Ir,x+w+dt,Ed0=ΨΔIl,x;Ir,x+d0,Edt=ΨΔIl,x+w;Ir,x+w+dt,where El and Er are the left and the right data term energy function of optical flow, respectively. Ed0 is the disparity term that represents the stereo matching at time t0. Edt is the disparity term at time t+1. Penalty function (Ψ(s2)=s2+ε2,ε=0.001) is introduced to guarantee convexity and differentiability of the whole data term energy function and remove the outliers at the same time. Outliers are the pixels that are not consistent with the pixels affected by noise, illumination change, and occlusion. γ can be determined according to the intensity of illumination change or just rely on the experiences.The most obvious thing should be noted here is that if there is an error in any term ofu,v,d0, and dt, the error would be propagated to others. We minimize the 4 data terms together to prevent the error from this situation.


## 3.3. Smooth Term Based on Anisotropic Total Variation

Huguet and Devernay [5] used flow-driven regularization that uses a robust and convex function to keep convexity and differentiability of the smooth term(7)Esmoothu,v,d0,dt=∫ΩΨ∇u2+∇v2+ρ∇d0-dt2+μ∇d02dxdy.Here Ψ is used to remove the outliers and preserve the discontinuities of the flow field. The magnitude of the gradient is(8)Ψ∇u=ux2+uy2+ε2.We further derive its scene flow form by minimizing the following problems. In order to simplify the equation, we only consider the 2D case(9)F=∫ΩΨ∇u2+∇v2dxdy=∫Ωux2+uy2+vx2+vy2+ε2dxdy.It corresponds to the diffusion equation as(10)∂F∂u=divΨ′∇u2+∇v2∇u=∂∂xuxux2+uy2+vx2+vy2+ε2+∂∂yuyux2+uy2+vx2+vy2+ε2,∂F∂v=divΨ′∇u2+∇v2∇v=∂∂xvxux2+uy2+vx2+vy2+ε2+∂∂yvyux2+uy2+vx2+vy2+ε2.In the directions of u and v, the weights are the same and are related to the sum of the squares of the gradient. We propose a new scheme that regularization includes information about the direction of the flow and we define it as anisotropic total variation smoothness term. In the 2D case,(11)TVanisotropic=∑i,j∇u+∑i,j∇v=∑i=1M∑j=1Nui,j-ui,j+1+ui,j-ui+1,j+∑i=1M∑j=1Nvi,j-vi,j+1+vi,j-vi+1,j.But for outliers, robustness function still should be used. Here we show that it is anisotropic. Similarly, in the 2D case, we define the following equation:(12)F=∫ΩΨux2+Ψuy2+Ψvx2+Ψvy2dxdy=∫Ωux2+ε2+uy2+ε2+vx2+ε2+vy2+ε2dxdy.We minimize it and obtain the following equation:(13)∂F∂u=∂∂xuxux2+ε2+∂∂yuyuy2+ε2,∂F∂v=∂∂xvxvx2+ε2+∂∂yvyvy2+ε2.Obviously, the weights only depend on the direction of each gradient of ux or vx and are not related to other direction, which can be regarded as anisotropic diffusion of flow field. In this way, we can control the rate of diffusion according to the gradient of flow field so that we can preserve the sharp edge at the position of nonflat area. And in the flat region, it has the same strength of smooth term. Then a new smooth term can be defined as(14)Esmoothu,v,d0,dt=∫ΩΨux2+Ψuy2+Ψvx2+Ψvy2+ρΨd0-dtx2+Ψd0-dty2+μΨd0x2+Ψd0y2dxdy,where ρ and μ are empirical parameters.


## 3.4. Smooth Term Based on Adaptive Anisotropic Total Variation Flow-Driven Method

A major deficiency of the total variation that Rudin et al. [17] pointed out is staircase in the steady state solution. Marquina and Osher [18] proposed a modified method in 2000, but they just relieved this phenomenon. Different from total variation, there is no such problem with L2 norm. Blomgren and Chan [13] used the following norm in the image restoration:(15)Rpu=∫Ω∇up∇udxdy,where p function monotonously decreases and meets the following requirement:(16)p∇u⟶2,∇u⟶01,∇u⟶∞.At the edge position, we have(17)Rpu≈∫Ω∇udxdy.In the flat region, we have(18)Rpu≈∫Ω∇u2dxdy.In other words, this norm can choose from L1 to L2 according to the gradient magnitude automatically. To achieve this, a rational form can be chosen as follows:(19)p∇u=∇u+2∇u+1.The Euler-Lagrange equation for (15) is(20)∂R∂u=divp∇u·∇up∇u-1+ln⁡p∇up′∇u·∇up∇u∇u∇u.Define a function to represent the equation of the bracket in (20):(21)q∇u=p∇u·∇up∇u-1+ln⁡p∇up′∇u·∇up∇u.Then (20) can be simplified as(22)∂R∂u=divq∇u∇u∇u.Equation (23) can be derived from (21):(23)q∇u=2∇u,∇u⟶01,∇u⟶∞.In the smooth region of flow, the regularization is similar to linear smoothing, and total variation model will be used closing to the flow edge. Combining the robust function Ψ(s2)=s2+ε2, our scene flow smooth term can be written as(24)Esmoothu,v,d0,dt=∫ΩΨux2pux+Ψuy2puy+Ψvx2pvx+Ψvy2pvy+ρΨd0-dtx2pd0-dtx+Ψd0-dty2pd0-dty+μΨd0x2pd0x+Ψd0y2pd0ydxdy.


## 3.5. Optimization

### 3.5.1. Euler-Lagrange Equations

Before the introduction of minimizing the energy function, we first introduce the following abbreviations for the left image. The similar abbreviation is also applicable to the right image.(25)Ilx≔∂xIlx+w,Ilxz≔∂xIlx+w-∂xIlx,Ily≔∂yIlx+w,Ilyz≔∂yIlx+w-∂yIlx,Ilz≔Ilx+w-Ilx,Ilyy≔∂yy2Ilx+w,Ilxx≔∂xx2Ilx+w,Ilxy≔∂xy2Ilx+w,Ilt+1≔Ilx+w.We also use abbreviations to replace the partial derivatives of the data term.(26)Ψ′fl≔∂xΨIlx+w-Ilx2+γ∇Ilx+w-∇Ilx2,Ψ′fr≔∂xΨIrx+w+dt-Irx+d02+γ∇Irx+w+dt-∇Irx+d02,Ψ′st≔∂xΨIrx+w+dt-Ilx+w2+γ∇Irx+w+dt-∇Ilx+w2.It should be noted that the regularization term we used is anisotropic, so the Euler-Lagrange equations are not the same as (21). Considering anisotropic total variation of individual variables(27)F=Ψux2pux+Ψuy2puy,we can get the following Euler equations:(28)∂F∂ux=2ux·pux·ψ′ux2·ψux2pux-1+ln⁡ψux2·p′ux·Ψux2puxuxux,∂F∂uy=2uy·puy·Ψ′uy2·Ψuy2puy-1+ln⁡Ψuy2·p′uy·Ψuy2puyuyuy.We define the q1 function as follows:(29)q1s=2s·ps·Ψ′s2·Ψs2ps-1+ln⁡Ψs2·p′s·Ψs2ps,and then (28) can be simplified as(30)∂F∂ux=q1uxuxux,∂F∂uy=q1uyuyuy.According to the variational principle, we minimize the energy function and achieve(31)∇Eu,v,d0,dt=0.Equation (31) corresponds to the 4 Euler equations. We calculate the partial derivation of E with respect to u and let it be equal to 0 to get the following equation:(32)βlΨ′fl·IlxIlz+γIlxxIlxz+IlxyIlyz+βrΨ′fr·IrxIrz+γIrxxIrxz+IrxyIryz+βstΨ′st·Irt+1-Ilt+1Irx-Ilx+γIrx-IlxIrxx-Ilxx+Iry-IlyIrxy-Ilxy-α∂∂xq1uxuxux+∂∂yq1uyuyuy=0.The partial derivatives of E with respect to v,d0-dt, and d0 produce similar equations as given below:(33)βlΨ′fl·IlyIlz+γIlxyIlxz+IlyyIlyz+βrΨ′fr·IryIrz+γIrxyIrxz+IryyIryz+βstΨ′st·Irt+1-Ilt+1Iry-Ily+γIrx-IlxIrxy-Ilxy+Iry-IlyIryy-Ilyy-α∂∂xq1vxvxvx+∂∂yq1vyvyvy=0,βrΨ′fr·IrxIrz+γIrxxIrxz+IrxyIryz+βstΨ′st·IrxIrt-Ilt+γIrxxIrx-Ilx+IrxyIry-Ily-αρ∂∂xq1dt-d0xdt-d0xdt-d0x+αρ∂∂yq1dt-d0ydt-d0ydt-d0y=0,βrΨ′fr·IrxtIrz+γIrxxtIrxz+IrxytIryz+βstΨ′st·IrxtIrz+γIrxxtIrxt-Ilxt+IrxytIryt-Ilyt-αρ+μ∂∂xq1d0xd0xd0x+∂∂yq1d0yd0yd0y=0.


### 3.5.2. Multiresolution Strategy

In order to solve the problem of large displacements in the scene and avoid local minima, we introduce the multiresolution method. Referring to Brox et al. [19], the sampling factor is set to 0.9 to avoid a larger span of adjacent layers and local minima. At the same time, we also have a serious problem in the scene flow estimation. If we start to calculate from low resolution layer, there will be a false disparity which is passed to the high resolution layer. For this reason, we start to calculate scene flow from middle resolution layer. A schematic diagram that uses the pyramid hierarchy is given in Figure 2.Figure 2
Multiresolution strategy for scene flow computation. The left side pyramid is the image at timet. The right side pyramid is the image at time t+1.In each level of the pyramid, we use the warping method. For the two frames of theK level, I0k is obtained according to a certain sampling standard. Itk can be obtained by motion compensation based on the image Itk-1 in upper level. We only need to minimize the residual of each layer of the scene flow and the error generated by the large displacement can be avoided. The evaluation of scene flow requires minimizing partial differential equations (32)-(33). This is done by using Euler-Lagrange equations. We use a nested fixed point iteration scheme suggested by Brox et al. [19] and linearize equations in each layer of pyramid. The inner loop is used to linearize brightness constancy assumption and outer loop linearizes the robust penalty function Ψ. We use SOR (Super Relaxation Iteration) as an iterative method. The full algorithm for scene flow estimation is given in the algorithm scheme (see Algorithm 1).Algorithm 1:Scene flow estimation.
Ensure: The images at time t and t+1 are in the same level Require: a>b≥c≥0 for level = c to b  do Left flow (ul,vl) adaptive anisotropic flow-driven Right flow (ur,vr) adaptive anisotropic flow-driven end for stereo (d0,dt) for level = b to a  do scene flow (ul,vl,d0,dt) adaptive anisotropic total variation flow-driven end for


## 3.5.1. Euler-Lagrange Equations

Before the introduction of minimizing the energy function, we first introduce the following abbreviations for the left image. The similar abbreviation is also applicable to the right image.(25)Ilx≔∂xIlx+w,Ilxz≔∂xIlx+w-∂xIlx,Ily≔∂yIlx+w,Ilyz≔∂yIlx+w-∂yIlx,Ilz≔Ilx+w-Ilx,Ilyy≔∂yy2Ilx+w,Ilxx≔∂xx2Ilx+w,Ilxy≔∂xy2Ilx+w,Ilt+1≔Ilx+w.We also use abbreviations to replace the partial derivatives of the data term.(26)Ψ′fl≔∂xΨIlx+w-Ilx2+γ∇Ilx+w-∇Ilx2,Ψ′fr≔∂xΨIrx+w+dt-Irx+d02+γ∇Irx+w+dt-∇Irx+d02,Ψ′st≔∂xΨIrx+w+dt-Ilx+w2+γ∇Irx+w+dt-∇Ilx+w2.It should be noted that the regularization term we used is anisotropic, so the Euler-Lagrange equations are not the same as (21). Considering anisotropic total variation of individual variables(27)F=Ψux2pux+Ψuy2puy,we can get the following Euler equations:(28)∂F∂ux=2ux·pux·ψ′ux2·ψux2pux-1+ln⁡ψux2·p′ux·Ψux2puxuxux,∂F∂uy=2uy·puy·Ψ′uy2·Ψuy2puy-1+ln⁡Ψuy2·p′uy·Ψuy2puyuyuy.We define the q1 function as follows:(29)q1s=2s·ps·Ψ′s2·Ψs2ps-1+ln⁡Ψs2·p′s·Ψs2ps,and then (28) can be simplified as(30)∂F∂ux=q1uxuxux,∂F∂uy=q1uyuyuy.According to the variational principle, we minimize the energy function and achieve(31)∇Eu,v,d0,dt=0.Equation (31) corresponds to the 4 Euler equations. We calculate the partial derivation of E with respect to u and let it be equal to 0 to get the following equation:(32)βlΨ′fl·IlxIlz+γIlxxIlxz+IlxyIlyz+βrΨ′fr·IrxIrz+γIrxxIrxz+IrxyIryz+βstΨ′st·Irt+1-Ilt+1Irx-Ilx+γIrx-IlxIrxx-Ilxx+Iry-IlyIrxy-Ilxy-α∂∂xq1uxuxux+∂∂yq1uyuyuy=0.The partial derivatives of E with respect to v,d0-dt, and d0 produce similar equations as given below:(33)βlΨ′fl·IlyIlz+γIlxyIlxz+IlyyIlyz+βrΨ′fr·IryIrz+γIrxyIrxz+IryyIryz+βstΨ′st·Irt+1-Ilt+1Iry-Ily+γIrx-IlxIrxy-Ilxy+Iry-IlyIryy-Ilyy-α∂∂xq1vxvxvx+∂∂yq1vyvyvy=0,βrΨ′fr·IrxIrz+γIrxxIrxz+IrxyIryz+βstΨ′st·IrxIrt-Ilt+γIrxxIrx-Ilx+IrxyIry-Ily-αρ∂∂xq1dt-d0xdt-d0xdt-d0x+αρ∂∂yq1dt-d0ydt-d0ydt-d0y=0,βrΨ′fr·IrxtIrz+γIrxxtIrxz+IrxytIryz+βstΨ′st·IrxtIrz+γIrxxtIrxt-Ilxt+IrxytIryt-Ilyt-αρ+μ∂∂xq1d0xd0xd0x+∂∂yq1d0yd0yd0y=0.


## 3.5.2. Multiresolution Strategy

In order to solve the problem of large displacements in the scene and avoid local minima, we introduce the multiresolution method. Referring to Brox et al. [19], the sampling factor is set to 0.9 to avoid a larger span of adjacent layers and local minima. At the same time, we also have a serious problem in the scene flow estimation. If we start to calculate from low resolution layer, there will be a false disparity which is passed to the high resolution layer. For this reason, we start to calculate scene flow from middle resolution layer. A schematic diagram that uses the pyramid hierarchy is given in Figure 2.Figure 2
Multiresolution strategy for scene flow computation. The left side pyramid is the image at timet. The right side pyramid is the image at time t+1.In each level of the pyramid, we use the warping method. For the two frames of theK level, I0k is obtained according to a certain sampling standard. Itk can be obtained by motion compensation based on the image Itk-1 in upper level. We only need to minimize the residual of each layer of the scene flow and the error generated by the large displacement can be avoided. The evaluation of scene flow requires minimizing partial differential equations (32)-(33). This is done by using Euler-Lagrange equations. We use a nested fixed point iteration scheme suggested by Brox et al. [19] and linearize equations in each layer of pyramid. The inner loop is used to linearize brightness constancy assumption and outer loop linearizes the robust penalty function Ψ. We use SOR (Super Relaxation Iteration) as an iterative method. The full algorithm for scene flow estimation is given in the algorithm scheme (see Algorithm 1).Algorithm 1:Scene flow estimation.
Ensure: The images at time t and t+1 are in the same level Require: a>b≥c≥0 for level = c to b  do Left flow (ul,vl) adaptive anisotropic flow-driven Right flow (ur,vr) adaptive anisotropic flow-driven end for stereo (d0,dt) for level = b to a  do scene flow (ul,vl,d0,dt) adaptive anisotropic total variation flow-driven end for


## 4. Experiments and Results

In our method, scene flow can be represented by optical flow(u,v) and disparity (d0,dt). We used data sets from Middlebury website to evaluate the accuracy of scene flow estimation. The datasets are captured from 8 calibrated cameras in each scene. The images from cameras 2 and 6 are taken as the stereo pairs at time t, while the images from cameras 4 and 8 are taken as the corresponding stereo pairs at time t+1. We evaluate the accuracy by calculating the RMS error in terms of optical flow, disparity map at time t, and disparity map at t+1 as well as the average angular error (AAE).We compared our results with Huguet and Devernay [5] and Basha et al. [6] using Venus, Cones, and Teddy from Middlebury datasets as shown in Figures 3, 4, and 5. Table 1 summarizes the RMS errors for optical flow, disparity at time t, and disparity at time t+1 to measure the deviation with the ground truth. AAE errors are computed to measure the deviation with the standard flow field.Table 1
Comparison of scene flow on Middlebury data set.
RMSO.F. RMSDisp. t RMSDisp. t+1 AAE Cones Huguet and Devernay 1.10 2.11 5.24 0.69 Basha et al. 0.58 2.48 2.49 0.39 Ours 0.79 1.22 2.71 0.56 Teddy Huguet and Devernay 1.25 2.27 6.93 0.51 Basha et al. 0.57 2.83 2.86 1.01 Ours 0.59 1.20 4.06 0.34 Venus Huguet and Devernay 0.31 0.97 1.48 0.98 Basha et al. 0.16 1.06 1.06 1.58 Ours 0.30 1.00 2.20 1.15Figure 3
Venus of Middlebury stereo datasets. (a) and (b) are left and right image in timet. (c) and (d) are ground truth of disparity maps. (e) and (f) are left and right image in time t+1. (g) is ground truth of optical flow in horizontal direction. (h) is ground truth of optical flow in vertical direction.
(a)
L 0 (b)
R 0 (c)
D 0 (d)
D t (e)
L 1 (f)
R 1 (g)
U (h)
VFigure 4
Cones of Middlebury stereo datasets. (a) and (b) are left and right image in timet. (c) and (d) are ground truth disparity maps. (e) and (f) are left and right image in time t+1. (g) is ground truth of optical flow in horizontal direction. (h) is ground truth of optical flow in vertical direction.
(a)
L 0 (b)
R 0 (c)
D 0 (d)
D t (e)
L 1 (f)
R 1 (g)
U (h)
VFigure 5
Teddy of Middlebury stereo datasets. (a) and (b) are left and right image in timet. (c) and (d) are ground truth disparity maps. (e) and (f) are left and right image in time t+1. (g) is ground truth of optical flow in horizontal direction. (h) is ground truth of optical flow in vertical direction.
(a)
L 0 (b)
R 0 (c)
D 0 (d)
D t (e)
L 1 (f)
R 1 (g)
U (h)
VRMS errors of optical flow (O.F.), disparity at timet (Disp. t), disparity at time t+1 (Dip. t+1), and AAE errors of optical flow are computed to evaluate the accuracy of our method.Table1 shows that the results of our method are more accurate than those of Huguet and Devernay [5] for RMS of optical flow and disparity. Comparing with Basha et al. [6], AAE errors of our method are more accurate. For the RMS errors in disparity t, our method is better than Huguet and Devernay and Basha et al.’s method. This benefits from proposed smooth term based on adaptive anisotropic total variational method.Figure6 shows the performance on the Venus dataset quantitatively. In terms of the component U, our method keeps more smooth and sharp boundary when comparing with Huguet and Devernay [5].Figure 6
Comparison of our U component and Huguet and Devernay’s [5]. (a) is the U component of Huguet and Devernay. (b) is our U component.
(a)
Huguet and Devernay [5] (b)
OursCones, Teddy, and Venus datasets have complex edges in the scenes. Traditional scene flow methods are suffering from the isotropic diffusion and receive lower flow accuracy. The proposed adaptive anisotropic diffusion smooth term can smooth scene flow along the direction of motion edges and chooseL1 or L2 regularization automatically, which can receive high flow accuracy at the position of edges and preserve motion edges.In order to compare our method with other scene flow methods, we used a rotating sphere dataset which includes relatively complex motion and consists of the hemispheres in opposite directions.Figure7 shows the left and right views of the sphere. The 3D motion u,v and disparity are also shown for comparison. Figure 8 shows results computed by our method. Table 2 summarizes the RMS errors of sphere sequence.Table 2
2D errors for rotating sphere sequence.
Huguet and Devernay [5] Wedel et al. [15] Vogel et al. [16] Ours RMS O.F. 0.69 0.77 0.63 0.63 RMS Disp. 3.80 10.90 2.84 3.67 The RMS errors of Huguet and Devernay [5], Wedel et al. [15], Vogel et al. [16], and our method are compared using optical flow (O.F.) and disparity map (Disp.).Figure 7
Rotating sphere which consists of the hemispheres in opposite directions. (a) and (b) are left and right images in time. (e) and (f) are left and right images in time. (c), (d), (g), and (h) are ground truth ofd0, dt, u, and v, respectively.
(a)
L 0 (b)
R 0 (c)
D 0 (d)
D t (e)
L 1 (f)
R 1 (g)
U (h)
VFigure 8
Rotating sphere which consists of the hemispheres in opposite directions. (a) and (b) ared0, dt computed by our method. (c) and (d) are u, v computed by our method.
(a)
L 0 (b)
R 0 (c)
D 0 (d)
D tTable2 shows that RMS error of our method in optical flow is the same as that of Vogel et al. [16] and better than others. For the RMS error in disparity at time t, our method produces a preferable result.Rotating sphere dataset has opposite directions motion, which is a very challenging task. The proposed method can smooth flow field along the motion edge and can also reduce the rate of diffusion at the position of motion edges. In this way, our method can preserve motion edges better.


## 5. Conclusion

In this article, we proposed an adaptive anisotropic total variation flow-driven regularization method in the variational framework for scene flow estimation. According to the corresponding relationship of the stereo image sequence, data term is divided into 4 parts and meets the brightness and gradient constancy assumption. Regularization term has the advantages of anisotropic smoothing, which applies different weights in the different directions.L1 and L2 norm can be chosen automatically according to the gradient of the flow field to preserve the boundary of scene flow and reduce the staircase effect at the same time. By minimizing the energy function, we obtain dense optical flow and disparity map. Pyramid optimization method and warping technology are also referenced to deal with large displacement problems.In the section of experiments, RMS errors of our method for Middlebury datasets Cones and Teddy are better than those of Huguet and Devernay [5] and close to those of Basha et al. [6]. For the rotation sphere dataset, our method’s result is better than those of Huguet and Devernay [5] and Wedel et al. [15] and has the same accuracy as that of Vogel et al. [16]. The experimental results show the effectiveness of our method.



---

*Source: 4908273-2018-05-20.xml*