Basha MVSF (multi-view scene flow)

Tali Basha (now Dekel) published a Multi-View Scene Flow (MVSF) at CVPR 2010, and IJCV 2013 [home] (with code) [author] [tr]. We search for possibilites for user interaction.

Algorithm

Given n cameras Ci with frames Ii at t=0 and t=1, and an initial estimate of Z0 in camera C0, the objective is to dermine the true depth Z0 and the 3D-motion U,V,W. A coarse-to-fine pyramid with η=.9 is used to remove the nonlinearity in the data term, an outer loop to remove the nonlinearity from the perspective division, and an inner loop to remove the nonlinearity from the nonquadratic penalizer ψ(s2)=[s22].5. In the algorithm sketch below, only U is used for illustration, V/W/Z0 go analoguously:

U0     for pyramid level "L"=14..0: unknowns helpers images maps coordinates
U1     U=0 or U=UL+1 U   I_0t, I_0tXY MapX_Id
U1     for outer iteration "O"=0..3:
U2         keep U fixed; rewarp Ii,0 and Ii,1 with current U; dU=0          dU, UXY   occl0t I_0tw, I_0twXY   MapX_I_0t   X, Xt, Xxyz, Xtxy, projxy0t,xyz, projxyt,uvw
U2         for inner iteration "I"=0..10:
U3             keep dU fixed; build Mx=B with current dU; x=0 dUXY occl0t, div0wz MapX_I_0tin Zin, Xin, Xtin
U3             for SOR iteration "S"=0..40:
U4                 keep ψ{flow,stereo1,stereo2} fixed; solve Mx=B x
U3             dU = x
U2         U += dU
U1     UL-1 = upsample(U)

Variables

Due to the 3D to 2D projection nature of the computation, a lot of coordinates are juggled around, as are their derivatives. The list below shows the variables in the MVSF algorithm in a sequenced overview. All equations refer to the technical report [tr] which has the most details. There, BC{m,s1,s2} → BC{flow,st.init,stereo} refers to brightness constancy and S{m,s} → S{flow,stereo} to smoothness.

U1     MapX_Id = Id(Ny,Nx)
U2                X = ComputeXY(Z0, MapX_Id) (eq.2)
         Xt = X+U      , Zt = Z0+W (eq.3,4)
         MapX_I_0 = ComputeProjMaps(X, Z0) (eq.5,16-17)
         MapX_I_t = ComputeProjMaps(Xt, Zt)
         occl0 = ComputeOcclusion(X, Z0, MapX_I_0) (eq.7, alg.1)
         occlt = ComputeOcclusion(Xt, Zt, MapX_I_t)     (eq.7)
         I_0w = Smooth(Remap(I_0, MapX_I_0)) (eq.15)
         I_tw = Smooth(Remap(I_t, MapX_I_t))
         I_0wX = Derivative(I_0w)
         I_twX = Derivative(I_tw)
         Z0X = GradientFast(Z0)
         UX = GradientFast(U)
         Xxyz = Gradient3DPoints(Z0, Z0X, MapX_Id)  , Xtxy = Xxy+UXY
         projx0,xyz = GradientProjCoor(X, Z0, Xxyz, Z0X)
         projxt,xyz = GradientProjCoor(Xt, Zt, Xtxyz, ZtX)
         projxt,uvw = GradientProjCoor(Xt, Zt, 01, 01)
U3              Zin = Z0+dZ0  , Ztin = Zin+W
             Xin = ComputeXY(Zin, MapX_Id)  , Xtin = Xin+U+dU
             MapX_I_0in = ComputeProjMaps(Xin, Zin)
             MapX_I_tin = ComputeProjMaps(Xtin, Ztin)
             dZ0X = GradientFast(dZ0)
             dUX = GradientFast(dU)
             div0wz = divergence(Z0X, dZ0X, UX, dUX) (eq.12,13)
             M, B, occl0t = MB(Z0, dZ0, U, dU, div0wz, occl0t, MapX_I_0tin, I_0(X), I_0tw(X), projx0t,xy, projxt,uvw)
U4                  x = SOR(M, B)

Mx=B

This is the system of linear equations to be solved and the best chance to insert user input into the optimization. Green entries mark optional parts (when solving for UVWZ instead of UVW). Note that the full M is actually a sparse (Nx*Ny*|uvwz|)2 matrix, and the M here represents the diagonal and neighbour locations. In the following, one 4-part for UVWZ is shown. The SOR solver in the section after this shows each of the UVWZ uses one row of M: U will use line0, V line1, W line2, Z line3.

M x = B
au1 av1 aw1 az1 bu bu bu bu dU bu
au2 av2 aw2 az2 bv bv bv bv . dV = bv
au3 av3 aw3 az3 bw bw bw bw dW bw
au4 av4 aw4 az4 bz bz bz bz dZ bz

where:

M0..3
au1 = α*μu * [ ∑i=1..4(bfi) ]  , bf2 = .5 * (div0[j,i] + div0[j,i-1]) , where bf1,bf2,bf3,bf4 are right(1),left(2),down(3),up(4) (Ss or Sm, eq.12,13)
+ ψflow * [ itwU ]2  , ψflow = (12 + [it+(itwU*dU)+(itwV*dV)+(itwW*dW) - i0 + (itwZ-i0wZ)*dZ]2)-.5 (BCm, Δti=1+, eq.11)
+ ψstereo * [ itwU - iltwU ]2  , ψstereo = (12 + [it+(itwU*dU)+(itwV*dV)+(itwW*dW)+(itwZ*dZ) -  ilt-(iltwU*dU)-(iltwV*dV)-(iltwW*dW)-(iltwZ*dZ)]2)-.5 (BCs2, Δ^i, eq.11)
+ ψflowL * [ iltwU ]2  , ψflowL = (12 + [ilt+(iltwU*dU)+(iltwV*dV)+(iltwW*dW)+(iltwZ*dZ) - il0]2)-.5 (BCm, Δti=0, eq.11)
aw2 = ψflow * [ itwW * itwV ]  , itwU = (itX*projxt,u)+(itY*projyt,u) , itX=(itwX*projyt,y - itwY*projyt,x)/det , det=projxt,x*projyt,y - projxt,y*projyt,x (J, eq.22)
+ ψstereo * [ (itwW - iltwW) * (itwV - iltwV) ]  , iltwU = (iltX*projxlt,u)+(iltY*projylt,u) , iltX=(iltwX*projylt,y - iltwY*projylt,x)/det , det=projxlt,x*projylt,y - projxlt,y*projylt,x (J, eq.22)
+ ψflowL * [ iltwW * iltwV ]  , it = Images_tw[1+], ilt=Images_tw[0] , i0=Images_0w[1+], il0=Images_0[0] , itwX=Images_tw_X[1+], iltwX=Images_tw_X[0]
az4 = αzz * [ ∑i=1..4(bfzi) ]  , i0wZ = (i0X*projx0,z)+(i0Y*projy0,z) (eq.18) , i0X=(i0wX*projy0,y - i0wY*projy0,x)/det , det=projx0,x*projy0,y - projx0,y*projy0,x (J, eq.22)
+ ψflow * [ itwZ - i0wZ ]2  , itwZ = (itX*projxt,z)+(itY*projyt,z) (eq.20) , itX=(itwX*projyt,y - itwY*projyt,x)/det , det=projxt,x*projyt,y - projxt,y*projyt,x (J, eq.22)
+ ψst.init * [ i0wZ ]2  , ψst.init = (12 + [i0+(i0wZ*dZ) - il0]2)-.5 (BCs1, Δi=1+, eq.11)
+ ψstereo * [ itwZ - iltwZ ]2
M4..7
bu = -α*μu * bf4  ,  bf4 = .5 * (div0[j,i] + div0[j-1,i]) , div0=(.0012 + [(UX+dUX)2+(UY+dUY)2+(VX+dVX)2+(VY+dVY)2])-.5 (Sm(uv), eq.13)
bw = ww * bfw2  , bfw2 = .5 * (divw[j,i] + divw[j,i-1]) , divw=(.0012 + [μw * ((WX+dWX)2+(WY+dWY)2)])-.5 (Sm(w), eq.13)
bz = zz * bfz1  , bfz1 = .5 * (divz[j,i] + divz[j,i+1]) , divz=(.0012 + [μz * ((Z0X+dZ0X)2+(Z0Y+dZ0Y)2)])-.5 (Ss, eq.12)
B
bu = -α*μu    * [ ∑i=1..4(bfi*bui) - ∑i=1..4(bfi) * bu0 ]  , bu0 = U[j,i], bu1=U[j,i+1], ...
- ψflow * [ itwU * (it-i0) ]
- ψstereo * [ (it-ilt) * (itwU-iltwU) ]
- ψflowL * [ iltwU * (ilt-il0) ]
bw = ww * [ ∑i=1..4(bfwi*bwi) - ∑i=1..4(bfwi) * bw0 ] , bw0 = W[j,i], bw1=W[j,i+1], bw2=W[j,i-1], ...
- ψflow * [ itwW * (it-i0) ]
- ψstereo * [ (it-ilt) * (itwW-iltwW) ]
- ψflowL * [ iltwW * (ilt-il0) ]
bz
=
zz
* [ ∑i=1..4(bfzi*bzi) - ∑i=1..4(bfzi) * bz0 ]  ,
bz0 =
Z0[j,i], ..., bz4=Z0[j-1,i]



-
ψflow * [ (itwZ-i0wZ) * (i0-it) ]






- ψst.init * [ (i0-il0) * i0wZ ]
- ψstereo * [ (it-ilt) * (itwZ-iltwZ) ]
- ψflowL * [ iltwZ * (ilt-il0) ]

Neighbours: The order 1..4 is always right(1),left(2),down(3),up(4), and the variables bf[i] and bu[i] always run concurrent, e.g. where bf2 asks left, bu2 asks left too. In M4..8, that order is reversed.

Observations:

SOR

Successive over-relaxation is a faster variant of the Jacobi iteration, using the diagonal and upper+lower triangular submatrixes. The row index i ∈ Nx*Ny*|uvwz| has separate entries for uvwz.

rnew = (1.-ΩSOR)*rold + ΩSOR * (b - sum)/diag = (1.-ΩSOR)*rold + ΩSOR * (bu - sumU)/au1
diag = Mi,0..3 where uvwz=0123 = au1, av2, aw3, or az4
b = Bi = bu, bv, bw, or bz
sumU = Mi,1*SORi+1 + Mi,2*SORi+2 + Mi,3*SORi+3 = av1*SORv + aw1*SORw + az1*SORz
+ Mi,4*SORi,1↑ + Mi,5*SORi,1 + Mi,6*SORi,1 + Mi,7*SORi,1 + bu*SORu,1↑ + bu*SORu,1 + bu*SORu,1 + bu*SORu,1
sumZ = Mi,0*SORi-3 + Mi,1*SORi-2 + Mi,2*SORi-1 = au4*SORu + av4*SORv + aw4*SORw
+ Mi,4*SORi,1↑ + Mi,5*SORi,1 + Mi,6*SORi,1 + Mi,7*SORi,1 + bz*SORz,1↑ + bz*SORz,1 + bz*SORz,1 + bz*SORz,1

For Jacobi iterations, set ΩSOR=1 so that rnew is not dependent on rold.

EOF (Sep:2014)