Numerical Approximation

The preceding Euler-Lagrange equations are nonlinear in their argument $\mathbf {w}=(u,v,1)^\top$. A first step towards a linear system of equations, which can be solved with common numerical methods, is the use of fixed point iterations on $\mathbf {w}$. In order to implement a multiscale approach, necessary to better approximate the global optimum of the energy, these fixed point iterations are combined with a downsampling strategy. Instead of the standard downsampling factor of $0.5$ on each level, it is proposed here to use an arbitrary factor $\eta\in(0,1)$, what allows smoother transitions from one scale to the next. Moreover, the full pyramid of images is used, starting with the smallest possible image at the coarsest grid. Let $\mathbf {w}^k = (u^k,v^k,1)^\top$, $k = 0,1,\ldots$, with the initialisation $\mathbf {w}^0=(0,0,1)^\top$ at the coarsest grid. Further, let $I^k_*$ be the abbreviations defined in (8) but with the iteration variable $\mathbf {w}^k$ instead of $\mathbf {w}$. Then $\mathbf {w}^{k+1}$ will be the solution of

$$\begin{array}{l} \Psi'((I^{k+1}_z)^2+\gamma((I^{k+1}_{xz})^2... ...ert^2)\nabla_{\hspace*{-0.5mm}3}v^{k+1}\right) =0. \end{array}$$ (9)

As soon as a fixed point in $\mathbf {w}^k$ is reached, we change to the next finer scale and use this solution as initialisation for the fixed point iteration on this scale.
Notice that we have a fully implicit scheme for the smoothness term and a semi-implicit scheme for the data term. Implicit schemes are used to yield higher stability and faster convergence. However, this new system is still nonlinear because of the nonlinear function $\Psi'$ and the symbols $I^{k+1}_*$. In order to remove the nonlinearity in $I^{k+1}_*$, first order Taylor expansions are used:

$$\begin{array}{lcl} I_z^{k+1} & \approx & I^k_z + I_x^k du^k ... ...& I_{yz}^k + I_{xy}^k du^k + I_{yy}^k dv^k\mbox{,} \end{array}$$

where $u^{k+1} = u^k + du^k$ and $v^{k+1} = v^k + dv^k$. So we split the unknowns $u^{k+1}$, $v^{k+1}$ in the solutions of the previous iteration step $u^k, v^k$ and unknown increments $du^k, dv^k$. For better readability let

$$\begin{array}{lcl} (\Psi')^k_{Data} &:=& \Psi'\Big((I^k_z + ... ...{\hspace*{-0.5mm}3} (v^{k}+dv^k)\vert^2) \mbox{,} \end{array}$$ (10)

where $(\Psi')^k_{Data}$ can be interpreted as a robustness factor in the data term, and $(\Psi')^k_{Smooth}$ as a diffusivity in the smoothness term. With this the first equation in system (9) can be written as

$$= (\Psi')^k_{Data}\cdot\Big(I^k_x \left(I^k_z + I_x^k du^k + I_y^k dv^k \right)\Big)$$ 0

$$+ \gamma\;(\Psi')^k_{Data}\cdot\Big(I^k_{xx} (I^k_{xz}+I_{xx}^k du^k+I_{xy}^k dv^k)+I^k_{xy}(I^k_{yz}+I_{xy}^k du^k+I_{yy}^k dv^k)\Big)$$

$$- \alpha\; \left((\Psi')^k_{Smooth}\nabla_{\hspace*{-0.5mm}3}(u^{k} + du^k)\right),$$ (11)

and the second equation can be expressed in a similar way. This is still a nonlinear system of equations for a fixed $k$, but now in the unknown increments $du^k, dv^k$. As the only remaining nonlinearity is due to $\Psi'$, and $\Psi$ has been chosen to be a convex function, the remaining optimisation problem is a convex problem, i.e. there exists a unique minimum solution.
In order to remove the remaining nonlinearity in $\Psi'$, a second, inner, fixed point iteration loop is applied. Let $du^{k,0}:=0$, $dv^{k,0}:=0$ be our initialisation and let $du^{k,l}, dv^{k,l}$ denote the iteration variables at some step $l$. Furthermore, let $(\Psi')^{k,l}_{Data}$ and $(\Psi')^{k,l}_{Smooth}$ denote the robustness factor and the diffusivity defined in (10) at iteration $k$, $l$. Then finally the linear system of equations in $du^{k,l+1}, dv^{k,l+1}$ reads

$$= (\Psi')^{k,l}_{Data}\cdot\Big(I^k_x \left(I^k_z + I_x^k du^{k,l+1} + I_y^k dv^{k,l+1} \right)$$ 0

$$+ \gamma I^k_{xx} (I^k_{xz}+I_{xx}^k du^{k,l+1} +I_{xy}^k dv^{k,l+1}) +\gamma I^k_{xy}(I^k_{yz}+I_{xy}^k du^{k,l+1} +I_{yy}^k dv^{k,l+1})\Big)$$

$$- \alpha\; \left((\Psi')^{k,l}_{Smooth}\nabla_{\hspace*{-0.5mm}3}(u^{k} + du^{k,l+1}) \right)$$ (12)

for the first equation. Using standard discretisations for the derivatives, the resulting sparse linear system of equations can now be solved with common numerical methods, such as Gauss-Seidel or SOR iterations. Expressions of type $I(\mathbf x + \mathbf {w^k})$ are computed by means of bilinear interpolation.