Differential Scene Flow from Light Field Gradients

Abstract

This paper presents novel techniques for recovering 3D dense scene flow, based on differential analysis of 4D light fields. The key enabling result is a per-ray linear equation, called the ray flow equation, that relates 3D scene flow to 4D light field gradients. The ray flow equation is invariant to 3D scene structure and applicable to a general class of scenes, but is under-constrained (3 unknowns per equation). Thus, additional constraints must be imposed to recover motion. We develop two families of scene flow algorithms by leveraging the structural similarity between ray flow and optical flow equations: local ‘Lucas–Kanade’ ray flow and global ‘Horn–Schunck’ ray flow, inspired by corresponding optical flow methods. We also develop a combined local–global method by utilizing the correspondence structure in the light fields. We demonstrate high precision 3D scene flow recovery for a wide range of scenarios, including rotation and non-rigid motion. We analyze the theoretical and practical performance limits of the proposed techniques via the light field structure tensor, a $3 \times 3$ matrix that encodes the local structure of light fields. We envision that the proposed analysis and algorithms will lead to design of future light-field cameras that are optimized for motion sensing, in addition to depth sensing.

Appendices

Appendix A: Proof of Result 2

Result 2 (Rank of structure tensor) Structure tensor $\mathbf {S}$ has three possible ranks: 0, 2, and 3 for a local 4D light field window. These correspond to scene patches with no texture (smooth regions), an edge, and 2D texture, respectively.

Proof

We first show the three cases of rank 0, 2 and 3, and then we prove that the structure tensor cannot be rank 1.

Since $rank(\mathbf {S})=rank(\mathbf {A}^T\mathbf {A})=rank(\mathbf {A})$, we only need to look at the rank of the $n\times 3$ matrix $\mathbf {A}$:

$$\begin{aligned} \mathbf {A}=\begin{bmatrix}\mathbf {A}_1&\mathbf {A}_2&\mathbf {A}_3\end{bmatrix} =\begin{bmatrix} L_{X1}&\quad L_{Y1}&\quad L_{Z1} \\ L_{X2}&\quad L_{Y2}&\quad L_{Z2} \\ \vdots&\quad \vdots&\quad \vdots \\ L_{Xn}&\quad L_{Yn}&\quad L_{Zn} \\ \end{bmatrix}. \end{aligned}$$

(16)

Case 1: Smooth region. In this case, $L_X = L_Y = L_Z = 0$ for all the locations in the light field window. Therefore, all the entries of $\mathbf {A}$ are zero, resulting in a rank 0 structure tensor. All three eigenvalues $\lambda _1, \lambda _2, \lambda _3 = 0$. As a result, it has a 3-D null space, and no motion vector can be recovered reliably.

Case 2: Single step edge. Without loss of generality, suppose the light field window corresponds to a fronto-parallel scene patch with a vertical edge, i.e., $L_Y = 0$ everywhere and thus $\mathbf {A}_2=\mathbf {0}$.

Consider a point P on the edge. Consider two rays from P that are captured on two horizontally-separated sub-aperture images indexed by $(x_a,y)$ and $(x_b,y)$. Let the coordinates of the two rays be $(x_a,y,u_a,v_a)$ and $(x_b,y,u_b,v_b)$, and let the light field gradients at these two rays be $(L_{Xa},L_{Ya},L_{Za})$ and $(L_{Xb},L_{Yb},L_{Zb})$. Note that these two gradients are two rows in $\mathbf {A}$. Recall that

$$\begin{aligned} L_{Z}=-\frac{u}{\varGamma }L_{X}-\frac{v}{\varGamma }L_{Y} . \end{aligned}$$

(17)

Since $L_{Ya}=0$ and $L_{Yb}=0$, we have,

$$\begin{aligned} L_{Za}=-\frac{u_a}{\varGamma }L_{Xa},\; L_{Zb}=-\frac{u_b}{\varGamma }L_{Xb}. \end{aligned}$$

(18)

Next, suppose there exists $k\ne 0$ such that $\mathbf {A}_3=k\mathbf {A}_1$. This implies

$$\begin{aligned} kL_{Xa}=-\frac{u_a}{\varGamma }L_{Xa},\;kL_{Xb}=-\frac{u_b}{\varGamma }L_{Xb}. \end{aligned}$$

(19)

By eliminating k from Eq. 19, we get $u_a=u_b$. However, since the scene point has a finite depth, the disparity $u_a - u_b \ne 0$. Therefore, $u_a \ne u_b$. This contradiction means such k doesn’t exist and $\mathbf {A}_1$ and $\mathbf {A}_3$ are linearly independent, which means the rank of $\mathbf {A}$ (and $\mathbf {S}$) is 2. As a result $\mathbf {S}$ has a 1-D null space (only one eigenvalue $\lambda _3 = 0$) and a 2D family of motions (motion orthogonal to the edge) can be recovered.

Case 3: 2D texture. In general $\mathbf {A}_1$, $\mathbf {A}_2$ and $\mathbf {A}_3$ are nonzero and independent. The structure tensor is full rank (rank $=3$) and the entire space of 3D motions are recoverable.

Now we show that the rank cannot be 1.

(Proof by contradiction) Assume there exists a 4D patch such that its corresponding matrix $\mathbf {A}$ is rank 1.

First $\mathbf {A}_1$ and $\mathbf {A}_2$ cannot both be zero. If we assume they are both zero then according to Eq. 17 all entries in $\mathbf {A}_3$ will also be zero, which results in a rank 0 matrix. Therefore $\mathbf {A}_1\ne \mathbf {0}$ or $\mathbf {A}_2\ne \mathbf {0}$.

Without loss of generality, assume $\mathbf {A}_1 \ne 0$. Since $\mathbf {A}$ is rank 1, there exists $k,l\in \mathbb {R}$ such that

$$\begin{aligned} \mathbf {A}_2= & {} k\mathbf {A}_1, \end{aligned}$$

(20)

$$\begin{aligned} \mathbf {A}_3= & {} l\mathbf {A}_1. \end{aligned}$$

(21)

Let us pick a ray $\mathbf {x_a}=(x_a,y_a,u_a,v_a)$ with light field gradient $(L_{Xa},L_{Ya},L_{Za})$ such that $L_{Xa}\ne 0$. Such $\mathbf {x_a}$ exists because $\mathbf {A}_1\ne 0$. Note that this ray is captured by the sub-aperture image indexed by $(x_a,y_a)$. Assume the scene point corresponding to $\mathbf {x_a}$ is observed in another sub-aperture image $(x_b,y_b)$ with $y_a=y_b$, in other words a sub-aperture image that is at the same horizontal line as $(x_a,y_a)$. Denote the corresponding ray as $\mathbf {x_b}=(x_b,y_a,u_b,v_b)$ with light field gradient $(L_{Xb},L_{Yb},L_{Zb})$. $L_{Xb}$ is also nonzero.

From Eq. 20 we know that $L_{Ya}=kL_{Xa}$ and $L_{Yb}=kL_{Xb}$. According to Eq. 17 we have

$$\begin{aligned} L_{Za}= & {} -\frac{u_a}{\varGamma }L_{Xa}-\frac{v_a}{\varGamma }\cdot L_{Ya}=-\frac{u_a+kv_a}{\varGamma }L_{Xa}, \end{aligned}$$

(22)

$$\begin{aligned} L_{Zb}= & {} -\frac{u_b}{\varGamma }L_{Xb}-\frac{v_b}{\varGamma }\cdot L_{Yb}=-\frac{u_b+kv_b}{\varGamma }L_{Xb}. \end{aligned}$$

(23)

From Eq. 21 we know that $L_{Za}=l L_{Xa}$, $L_{Zb}=l L_{Xb}$, and from Eqs. 22–23 we have

$$\begin{aligned} l=-\frac{u_a+kv_a}{\varGamma }=-\frac{u_b+kv_b}{\varGamma }. \end{aligned}$$

(24)

However since $x_a\ne x_b$, $y_a=y_b$, we have $u_a\ne u_b$ and $v_a=v_b$ due to simple epipolar geometry. Therefore Eq. 24 cannot hold, which means our assumption is false, and $rank(\mathbf {A})$ cannot be 1. $\square$

Appendix B: Implementation Details

1.1 B.1 Global Method

In Sect. 5, we introduced the global ‘Horn–Schunck’ ray flow method, which solves for the 3D scene motion by minimizing a functional:

$$\begin{aligned} E(\mathbf {V})= & {} \underbrace{E_D (\mathbf {V})}_\mathbf{Error term } + \underbrace{E_S (\mathbf {V})}_\mathbf{Smoothness term } , \ \ \text {where}\nonumber \\ E_D (\mathbf {V})= & {} \int _{\varOmega }\left( L_X V_X + L_Y V_Y + L_Z V_Z + L_t \right) ^2 dx\,dy\,du\,dv , \nonumber \\ E_S (\mathbf {V})= & {} \int _{\varOmega } \left( \lambda |\nabla V_X|^2+\lambda |\nabla V_Y|^2+\lambda _Z|\nabla V_Z|^2 \right) dx\,dy\,du\,dv .\nonumber \\ \end{aligned}$$

(25)

This is a convex functional and its minimum can be found by the Euler–Lagrange equations,

$$\begin{aligned} L_X(L_XV_X+L_YV_Y+L_ZV_Z)-\lambda \varDelta V_X&=-L_XL_t, \nonumber \\ L_Y(L_XV_X+L_YV_Y+L_ZV_Z)-\lambda \varDelta V_Y&=-L_YL_t , \nonumber \\ L_Z(L_XV_X+L_YV_Y+L_ZV_Z)-\lambda _Z\varDelta V_Z&=-L_ZL_t. \end{aligned}$$

(26)

These equations can be discretized as a sparse linear system, and is solved using Successive Over-Relaxation (SOR).

1.2 B.2 Structure-Aware Global Method

In this section we discuss an enhanced version of the structure-aware global method, which adopts the enhancing techniques for the local and global method, as discussed in Sect. 6 of the main paper.

1.2.1 B.2.1 Data Term

The data term is defined as:

$$\begin{aligned} E_D (\mathbf {V})=\int _{\varOmega _c}\sum _{\mathbf {x}_i\in \mathscr {P}(u,v)}h_i\rho _D((L_0(\mathbf {x_i})-L_1(\mathbf {w}(\mathbf {x_i},\mathbf {V})))^2)du\,dv, \end{aligned}$$

(27)

where $\mathscr {P}(u,v)$ is the 2D plane defined in Equation (11) in the main paper.

Weighted 2D window Notice that each ray in the 2D plane is given a different weight $h_i$, which is given by

$$\begin{aligned}&h_i=h(\mathbf {x_i},\mathbf {x_c})=h_g(\mathbf {x_i},\mathbf {x_c})\cdot h_o(\mathbf {x_i},\mathbf {x_c}), \end{aligned}$$

(28)

$$\begin{aligned}&h_g(\mathbf {x_i},\mathbf {x_c}) = e^{-\frac{(x_i-x_c)^2+(y_i-y_c)^2+(u_i-u_c+\alpha (x_i-x_c))^2+(v_i-v_c+\alpha (y_i-y_c))^2}{\sigma _g^2}},\nonumber \\\end{aligned}$$

(29)

$$\begin{aligned}&h_o(\mathbf {x_i},\mathbf {x_c}) = e^{-\frac{(d_{\alpha i}-d_{\alpha c})^2}{\sigma _o^2}}, \end{aligned}$$

(30)

where $\mathbf {x_c}$ denotes the center ray of the window. $d_\alpha =1/\alpha$ and is proportional to the actual depth of the scene point.

$h_g$ defines a Gaussian weight function that is based on the distance between $\mathbf {x_i}$ and $\mathbf {x_c}$ in the 2D plane. $h_o$ defines an occlusion weight by penalizing the difference in the estimated disparity $\alpha$ at $\mathbf {x_i}$ and $\mathbf {x_c}$. Notice that not all rays on $\mathscr {P}(u,v)$ corresponds to the same scene point as $\mathbf {x}_c$ because of occlusion. If the scene point corresponding to $\mathbf {x}_i$ occludes or is occluded by the scene point corresponding to $\mathbf {x}_c$, they will have a different $\alpha$ and thus a small value of $h_o$.

1.2.2 B.2.2 Smoothness Term

The smoothness term is defined as

$$\begin{aligned} E_S(\mathbf {V})&=\int _{\varOmega }g(\mathbf {x})(\lambda \sum _{i=1}^2\rho _S(V_{X(i)}^2) +\lambda \sum _{i=1}^2\rho _S(V_{Y(i)}^2)\nonumber \\&\quad +\lambda _Z\sum _{i=1}^2\rho _S(V_{Z(i)}^2))du\, dv, \end{aligned}$$

(31)

where $V_{X(i)}$ is short for $\frac{\partial V_X}{\partial u^{(i)}}$. (For simplicity we denote u, v as $u^{(1)},u^{(2)}$ respectively.) $g(\mathbf {x})$ is a weight function that varies across the light field. The error term $E_C(\mathbf {V})$ uses the warp function (Eq. 8) in the main paper.

Best practices from optical flow We choose the penalty function $\rho _D$, $\rho _S$ to be the generalized Charbonnier penalty function $\rho (x^2)=(x^2+\epsilon ^2)^a$ with $a=0.45$ as suggested in Sun et al. (2010).

Weight function for the regularization term The weight function $g(\mathbf {x})$ consists of two parts, which is combined using a harmonic mean:

$$\begin{aligned} g(\mathbf {x}) = \frac{g_c(\mathbf {x})g_d(\mathbf {x})}{g_c(\mathbf {x})+g_d(\mathbf {x})}. \end{aligned}$$

(32)

Consistency between XY-motion and Z-motion. In practice we notice that motion discontinuity is preserved better in XY-motion than in Z-motion. To improve the accuracy of the Z-motion, we solve the 3D motion $\mathbf {V}$ in a two-step process. We compute an initial estimate of the XY-motion, denoted as $\mathbf {U}=(U_X,U_Y)$, in the first pass. We then use $\mathbf {U}$ to compute a weight map for the regularization term:

$$\begin{aligned} g_c(\mathbf {x}) = \frac{1}{1+(|\nabla U_X|^2+|\nabla U_Y|^2)/\sigma _c^2}, \end{aligned}$$

(33)

where $\mathscr {N}(\mathbf {x})$ denotes a local neighborhood of the point $\mathbf {x}$. Then the full 3D motion $\mathbf {V}$ is computed in a second pass. Notice that $g(\mathbf {x})$ is small where gradient of $\mathbf {U}$ is large, in other words the regularization term will contribute less to the whole energy where there is a discontinuity in $\mathbf {U}$.

Consistency between motion boundaries and depth boundaries. We also assume the motion boundaries are likely to align with depth boundaries. In other words, we give a lower weight for points where the depth gradient is large:

$$\begin{aligned} g_d(\mathbf {x}) = \frac{1}{1+|\nabla d_\alpha |^2/\sigma _d^2}. \end{aligned}$$

(34)

1.2.3 B.2.3 Optimization

The error term $E_D(V)$ can be linearized as ,

$$\begin{aligned} E_D' (\mathbf {V})=\int _{\varOmega _c}\sum _{\mathbf {x}_i\in \mathscr {P}(u,v)}h_i\rho _D((L_{Xi}V_X+L_{Yi}V_Y+L_{Zi}V_Z+L_{ti})^2)du\,dv. \end{aligned}$$

(35)

Then the entire energy $E'=E_D'+E_S$ can be minimized using Euler–Lagrange equations:

$$\begin{aligned}&\sum _{\mathbf {x}_i\in \mathscr {P}(u,v)}h_i\rho _D'L_X\delta _L-\lambda \sum _{i=1}^2\frac{\partial }{\partial u^{(i)}}(g\rho _S'(V_{X(i)})V_{X(i)})= \nonumber \\&\quad -\sum _{\mathbf {x}_i\in \mathscr {P}(u,v)}h_i\rho _D'L_XL_t, \nonumber \\&\sum _{\mathbf {x}_i\in \mathscr {P}(u,v)}h_i\rho _D'L_Y\delta _L-\lambda \sum _{i=1}^2\frac{\partial }{\partial u^{(i)}}(g\rho _S'(V_{Y(i)})V_{Y(i)})=\nonumber \\&\quad -\sum _{\mathbf {x}_i\in \mathscr {P}(u,v)}h_i\rho _D'L_YL_t, \nonumber \\&\sum _{\mathbf {x}_i\in \mathscr {P}(u,v)}h_i\rho _D'L_Z\delta _L-\lambda _Z\sum _{i=1}^2\frac{\partial }{\partial u^{(i)}}(g\rho _S'(V_{Z(i)})V_{Z(i)})=\nonumber \\&\quad -\sum _{\mathbf {x}_i\in \mathscr {P}(u,v)}h_i\rho _D'L_ZL_t , \end{aligned}$$

(36)

where $\rho _D'$ is short for $\rho _D'((L_XV_X+L_YV_Y+L_ZV_Z+L_t)^2)$, $\delta _L=L_XV_X+L_YV_Y+L_ZV_Z$. Again, these equations are discretized and solved using SOR. The linearization step can then be repeated in an iterative, multi-resolution framework.