Error Analysis and Condition Estimation of the Pyramidal Form of the Lucas-Kanade Method in Optical Flow

Abstract

Optical flow is the apparent motion of the brightness patterns in an image. The pyramidal form of the Lucas-Kanade ($\mathrm{LK}$) method is frequently used for its computation but experiments have shown that the method has deficiencies. Problems arise because of numerical issues in the least squares ($\mathrm{LS}$) problem $min{∥Ax-b∥}_2^2$, $A\in {\mathbb{R}}^{m\times 2}$ and $m\gg 2$, which must be solved many times. Numerical properties of the solution $x_0=A^{†}b$ = ${\left(A^{T}A\right)}^{-1}{A}^{T}b$ of the $\mathrm{LS}$ problem are considered and it is shown that the property $m\gg 2$ has implications for the error and stability of $x_0$. In particular, it can be assumed that b has components that lie in the column space (range) $\mathcal{R}\left(A\right)$ of A, and the space that is orthogonal to $\mathcal{R}\left(A\right)$, from which it follows that the upper bound of the condition number of $x_0$ is inversely proportional to $cos\theta$, where $\theta$ is the angle between b and its component that lies in $\mathcal{R}\left(A\right)$. It is shown that the maximum values of this condition number, other condition numbers and the errors in the solutions of the $\mathrm{LS}$ problems increase as the pyramid is descended from the top level (coarsest image) to the base (finest image), such that the optical flow computed at the base of the pyramid may be computationally unreliable. The extension of these results to the problem of total least squares is addressed by considering the stability of the optical flow vectors when there are errors in A and b. Examples of the computation of the optical flow demonstrate the theoretical results, and the implications of these results for extended forms of the $\mathrm{LK}$ method are discussed.

1. Introduction

Optical flow is the apparent motion of the brightness patterns in an image. This is different from the motion field, which is the projection onto two dimensions of motion in a three dimensional environment. Both these motions are usually termed optical flow even though they are not necessarily the same ([§2] [1]). There are several methods for the computation of optical flow and it is assumed in all the methods that it is locally smooth. The methods can be classified as global, for example, the Horn-Schunck method [2], or local, for example, the Lucas-Kanade ($\mathrm{LK}$) method [3], and a review of these methods, and other methods, is in [4,5,6,7,8]. Local methods are more robust to noise but they do not yield a dense optical flow field, and global methods are more sensitive to noise but they yield a dense optical flow field. The computation of the optical flow has been the focus of much research, particularly in images that are corrupted by significant noise, or images in which two or more objects move, or images in which there are temporal changes in illumination [9].

Comparisons of the different methods for the computation of optical flow are in [1,10]. The conclusions from these comparisons are empirical and based on quantitative measures of the errors in the computed results, but numerical and linear algebraic considerations of the methods are not included in these comparisons. These issues are addressed in this paper by error analysis and condition estimation of the pyramidal form of the $\mathrm{LK}$ method. This method requires that many least squares ($\mathrm{LS}$) problems be solved at each level of the pyramid and it is shown that, even in the absence of occlusions, motion discontinuities and other problematic phenomena, the maximum values of the errors and condition numbers of the solutions of the $\mathrm{LS}$ problems (the optical flow) increase as the pyramid is descended from the top level to the bottom level.

Kearney, Thompson and Boley [11] consider condition estimation and error analysis for the computation of the optical flow from two points with pixel coordinates $p_i=(x_i,y_i)$ and $p_j=(x_j,y_j)$ at time t, where $p_i$ is near $p_j$ such that a first order Taylor approximation is appropriate. This approximation yields

$$\begin{array}{c}\hfill G\omega =-b,\end{array}$$

(1)

where

$$\begin{array}{c}\hfill G=\left[\begin{array}{cc}{I}_x^{\left(i\right)}& I_y^{\left(i\right)}\\ I_x^{\left(j\right)}& I_y^{\left(j\right)}\end{array}\right],\omega =\left[\begin{array}{c}u\\ v\end{array}\right],b=\left[\begin{array}{c}{I}_t^{\left(i\right)}\\ I_t^{\left(j\right)}\end{array}\right],\end{array}$$

$$\begin{array}{c}\hfill I_x^{\left(i\right)}=\frac{\partial I^{\left(i\right)}}{\partial x},I_y^{\left(i\right)}=\frac{\partial I^{\left(i\right)}}{\partial y},I_t^{\left(i\right)}=\frac{\partial I^{\left(i\right)}}{\partial t},\end{array}$$

$\omega$ is the optical flow and $I^{\left(i\right)}=I^{\left(i\right)}(x,y,t)$ is the intensity of the ith pixel. The work in this paper follows the analysis in [11] because issues of computational linear algebra, specifically condition estimation and error analysis, are addressed. There are, however, differences because an arbitrary number of points, rather than two points, is considered in this paper since the $\mathrm{LK}$ method yields an overdetermined set of equations, rather than the problem (1), for which the coefficient matrix is square and of order two. It is therefore necessary to compute $x_0$, the optical flow,

$$\begin{array}{c}\hfill x_0=arg\underset{x}{min}{∥Ax-b∥}^2,A\in {\mathbb{R}}^{m\times n},\end{array}$$

(2)

where $∥·∥={∥·∥}_2$, $m>n$ and rank $A=n$, for each window at each level of the pyramid. The optical flow problem requires that $n=2$ but it is assumed in the theoretical analysis that $n\le m$ is arbitrary.

The first and second columns of A contain, respectively, the intensity gradients with respect to x and y, and b contains the negative of the temporal intensity gradients, and errors are therefore present in A and b, such that it is necessary to solve the total least squares ($\mathrm{TLS}$) problem ([§6.3] [12]),

$$\begin{array}{c}\hfill \begin{array}{cc}\left\{x_{tls},\delta A_{tls},\delta b_{tls}\right\}=arg\underset{x,\delta A,\delta b}{min}{∥\left[\begin{array}{cc}\delta A\hfill & \delta b\hfill \end{array}\right]∥}_F\hfill & \\ & \\ subjectto\left(A+\delta A\right)x=b+\delta b,\hfill \end{array}\end{array}$$

(3)

where ${∥·∥}_F$ denotes the Frobenius norm. The method of $\mathrm{TLS}$ has been used for the computation of optical flow [13], but an estimate of its numerical stability, which is a measure of its sensitivity to noise, has not been established. This measure is called the condition number and this paper addresses this issue by considering the condition number of the computed optical flow. This condition number must be derived for the solution $x_{tls}$ of the $\mathrm{TLS}$ problem (3), but it is shown that it is first necessary to consider the condition number of the solution $x_0$ of the $\mathrm{LS}$ problem (2).

The results of this paper are summarised:

• It is shown that the maximum value of the error in the solution $x_0$ of (2) increases as the pyramid is descended from the top level to the lowest level.
• The condition number $\kappa \left(A\right)=∥A∥∥A^{†}∥$ of A, where $A^{†}={\left(A^{T}A\right)}^{-1}{A}^T$, is a function of A only but the solution $x_0$ of the $\mathrm{LS}$ problem is a function of A and b. A refined condition number, called the effective condition number $\eta (A,b)$ that is a function of A and b, is therefore developed.
• An upper bound of $\eta (A,b)$ is $\raisebox{1ex}{$\kappa \left(A\right)$}\!\left/ \!\raisebox{-1ex}{$cos\theta $}\right.$, where $\theta$ is the angle between b and its component that lies in the column space (range) $\mathcal{R}\left(A\right)$ of A. It follows that even if A is well conditioned, the solution $x_0$ may be ill conditioned if $\theta$ is large.
• The dimensions of the coefficient matrix A satisfy $m\gg n$ and it can therefore be assumed that a significant component of b lies in the space that is orthogonal to $\mathcal{R}\left(A\right)$. It is shown that the error in $x_0$ is therefore large and that there is a simple relationship between this error and $cos\theta$. It is concluded that the angle $\theta$ arises in the error and condition numbers of $x_0$, and it is therefore an important measure to consider when assessing the computational reliability of the optical flow.
• It is shown that the maximum values of the condition numbers and errors in the solutions of the $\mathrm{LS}$ problems (one $\mathrm{LS}$ problem for each window at each level of the pyramid) increase as the pyramid is descended from the coarsest image to the finest image. This shows that the advantage of the pyramidal implementation of the $\mathrm{LK}$ method - the ability to cope with large displacements between successive images - must be balanced against the numerical measures that quantify the robustness and reliability of the solutions of the many $\mathrm{LS}$ problems that are solved in the $\mathrm{LK}$ method.
• The condition number $\raisebox{1ex}{$\kappa \left(A\right)$}\!\left/ \!\raisebox{-1ex}{$cos\theta $}\right.$ of $x_0$ is limited because it quantifies the effect on $x_0$ of errors in b, but errors in A are not considered. A condition number $\rho (A,b)$ of $x_0$ that includes errors in A and b is developed because it is more appropriate for the $\mathrm{LK}$ method. This is a non-linear condition number, in contrast to the linear condition number $\raisebox{1ex}{$\kappa \left(A\right)$}\!\left/ \!\raisebox{-1ex}{$cos\theta $}\right.$, and it is shown that the dominant term in the expression for $\rho (A,b)$ is $\kappa \left(A\right)\eta (A,b)$.

The method of least median of squares $\left(\mathrm{LMedS}\right)$ is used in [14] to solve simultaneously the problems of temporal variation of illumination and motion discontinuities induced by the relative motion of two or more objects. The $\mathrm{LMedS}$ method requires the solution of the $\mathrm{LS}$ problem (2), and the pyramidal form of the $\mathrm{LK}$ method is also used in feature tracking [15]. More generally, the theory and results in this paper are applicable to all problems in which the $\mathrm{LK}$ method is used. There are several forms of the $\mathrm{LK}$ method and they differ in the expressions used for the calculation of the derivatives of the intensity. For example, one form of the $\mathrm{LK}$ method requires that A and b in (2) are premultiplied by a diagonal matrix, such that a larger weight is assigned to pixels in the centre of each window, and a smaller weight is assigned to pixels at its borders ([§4.3.2] [16]).

The $\mathrm{LK}$ method is described in Section 2 and the error in the solution $x_0$ of the $\mathrm{LS}$ problem is considered in Section 3. It is shown that there is a simple relationship between $cos\theta$ and the error, and an expression for $cos\theta$ is developed in Section 4. Linear and non-linear condition numbers that consider errors in b, and errors in A and b, respectively, of the optical flow are considered in Section 5. Examples that show the errors and condition numbers of the optical flow at each level of the pyramid are in Section 6. The paper is summarised in Section 7 and it is shown that the numerical problems with the solution of the $\mathrm{LS}$ problem that are addressed in this paper provide the mathematical motivation for the generalised form of the $\mathrm{LK}$ method that includes warping functions [17].

2. The Lucas-Kanade Method

The $\mathrm{LK}$ method is derived from the assumption of constant intensity $I(x,y,t)$ as a pixel at position $(x,y)$ at time t moves to position $(x+\delta x,y+\delta y)$ at time $t+\delta t$, such that

$$\begin{array}{c}\hfill I(x,y,t)=I(x+\delta x,y+\delta y,t+\delta t),\end{array}$$

to first order. It therefore follows that

$$\begin{array}{c}\hfill \frac{\partial I}{\partial x}\delta x+\frac{\partial I}{\partial y}\delta y+\frac{\partial I}{\partial t}\delta t=0,\end{array}$$

(4)

and thus if

$$\begin{array}{c}\hfill u=\frac{\partial x}{\partial t}andv=\frac{\partial y}{\partial t},\end{array}$$

then

$$\begin{array}{c}\hfill I_{x}u+I_{y}v+I_t=0,\end{array}$$

(5)

where

$$\begin{array}{c}\hfill I_x=\frac{\partial I}{\partial x},I_y=\frac{\partial I}{\partial y}and{I}_t=\frac{\partial I}{\partial t}.\end{array}$$

Equation (5) is one equation in two unknowns (the components u and v of the optical flow), and it therefore possesses an infinite number of solutions. This problem of a non-unique solution is addressed by considering a square window in each image at times t and $t+\delta t$, and if each window contains p points, then the application of (5) to each of these points yields,

$$\begin{array}{c}\hfill I_x^{\left(i\right)}u+I_y^{\left(i\right)}v+I_t^{\left(i\right)}=0,i=1,\cdots ,p,\end{array}$$

which can be combined into one equation,

$$\begin{array}{c}\hfill Ax=b,A\in {\mathbb{R}}^{p\times 2},x\in {\mathbb{R}}^{2}andb\in {\mathbb{R}}^p.\end{array}$$

(6)

Equation (6) is solved in the $\mathrm{LS}$ sense (2), which yields the solution,

$$\begin{array}{c}\hfill x_0=A^{†}b,A^{†}={\left(A^{T}A\right)}^{-1}{A}^T,\end{array}$$

(7)

where $x_0={\left[uv\right]}^T$ is the optical flow. It follows from (6) that the motion in each window is a translation only, where u and v represent its horizontal and vertical components, respectively.

The solution $x_0$ assumes there are errors in b only, and that A does not contain errors. The qualifications that follow from this assumption must be addressed because A and b contain, respectively, the spatial and temporal derivatives of the intensity $I(x,y,t)$. The inclusion of errors in A and b requires that the $\mathrm{TLS}$ problem be considered ([§6.3] [12]), and the effect of perturbations in A and b on the optical flow $x_0$ is considered by the development of a non-linear condition number of $x_0$ in Section 5.2.

The computation of the translation $x_0$ between two images of a dynamic scene occurs in several applications, including remote sensing and medical imaging, and it has been shown that noise and the first order approximation (4) introduce bias in the estimate of $x_0$ [11,18,19,20]. Robinson and Milanfar [20] use prior knowledge of the image and the translation to design gradient-estimation filters but the bias due to noise is not considered, and thus the $\mathrm{LK}$ method yields poor results for images whose signal-to-noise ratio is less than about 25 dB. Pham et al. [19] reduce the bias by iteratively computing the translation and upsampling the matrices that store the components u and v of the optical flow, which is the pyramid procedure proposed by Lucas and Kanade [3], and refined by Baker and Matthews [17]. Pham et al. [19] report that a pyramid with three levels yields an optical flow field with a very small bias.

The first order approximation (4) requires that

$$\begin{array}{c}\hfill \left|\delta x\right|,\left|\delta y\right|,\left|\delta t\right|\ll 1,\end{array}$$

(8)

and if this condition is not satisfied, a pyramid is formed in which the images at the ith level are formed by passing the images at the $(i-1)$th level through a low pass filter and then downsampling the filtered images. The coefficients of the low pass filter are [21,22],

$$\begin{array}{c}\hfill \left(\frac{1}{16}\right)\left[\begin{array}{c}1\\ 4\\ 6\\ 4\\ 1\end{array}\right]\left(\frac{1}{16}\right)\left[\begin{array}{ccccc}1& 4& 6& 4& 1\end{array}\right],\end{array}$$

(9)

which is approximately equal to a Gaussian filter with standard deviation $\sigma =1$.

The pyramidal implementation is shown in Figure 1 for two images F and G that are taken at times t and $t+\delta t$ respectively. The optical flow is computed at the top of the pyramid, level 4 (the coarsest image) in Figure 1, which is the initial condition for the computation of the optical flow at level 3 of the pyramid. It follows from (8) that (5) yields good estimates of u and v if the translation between the images at times t and $t+\delta t$ is small. Only the increments between successive levels of the pyramid are therefore calculated, after warping the two images based on the flow field estimate from the coarse level. This operation yields two matrices, one for each component of the optical flow, and each matrix is upsampled to yield matrices of twice the number of rows and columns, and then passed through a modified form of the filter (9),

$$\begin{array}{c}\hfill \left(\frac{1}{8}\right)\left[\begin{array}{c}1\\ 4\\ 6\\ 4\\ 1\end{array}\right]\left(\frac{1}{8}\right)\left[\begin{array}{ccccc}1& 4& 6& 4& 1\end{array}\right].\end{array}$$

(10)

This process is continued to the base of the pyramid, level 1, and the optical flow at this level is the desired optical flow.

The difference in the scale factors for the filters for downsampling the images and upsampling the optical flow field arises because all the pixels in the images that are downsampled are included in the computation, but rows and columns of zeros are added to the matrices that store the components u and v before they are upsampled, such that a matrix of order $M\times N$ is increased to a matrix of order $2M\times 2N$. The computation of all the pixel values in these larger matrices cannot include the added rows and columns of zeros and thus only one quarter of the pixels in the matrices of order $2M\times 2N$ are included in the computation.

The matrices of the components of the optical flow must be upsampled by the filter (10) and then multiplied by two in order to preserve the values of its components between successive levels of the pyramid. In particular, two adjacent pixels have unit separation before upsampling, but their separation increases to two after upsampling. The preservation of the values of the components of the optical flow in the upsampled images therefore requires that they be multiplied by two.

The discussion above shows that the pyramidal form of the $\mathrm{LK}$ method requires that the $\mathrm{LS}$ problem be solved many times at each level of the pyramid, and thus the error in its solution must be considered. This issue is addressed in Section 3 and it is shown in Section 4 that this error is closely related to the value of the angle $\theta$ between b and its component that lies in the column space $\mathcal{R}\left(A\right)$ of A. The computed optical flow must be numerically stable and it is therefore necessary to consider its numerical condition. This issue is addressed in Section 5, and linear and non-linear condition numbers are considered.

3. The Error in the Solution of the $\mathrm{LS}$ Problem

This section considers the error e in the solution of the $\mathrm{LS}$ problem (7),

$$\begin{array}{c}\hfill e=\frac{∥b-A{x}_0∥}{∥b∥}=\frac{∥(I-A{A}^{†})b∥}{∥b∥}=\frac{∥(I-\Sigma {\Sigma }^{†})c∥}{∥c∥},\end{array}$$

(11)

where the singular value decomposition of A is $U\Sigma V^T$, U and V are orthogonal matrices of orders m and n respectively,

$$\begin{array}{c}\hfill \Sigma =\left[\begin{array}{c}{\Sigma }_1\\ 0\end{array}\right]\in {\mathbb{R}}^{m\times n},{\Sigma }_1=diag\left[\begin{array}{cccc}{\sigma }_1& {\sigma }_2& \cdots & {\sigma }_n\end{array}\right],\end{array}$$

the singular values ${\sigma }_i$ satisfy ${\sigma }_i\ge {\sigma }_{i+1}$, $i=1,\cdots ,n-1,$ and

$$\begin{array}{c}\hfill c=U^{T}b.\end{array}$$

(12)

The matrix U is partitioned into two matrices $U_1$ and $U_2$,

$$\begin{array}{c}\hfill U=\left[\begin{array}{cc}{U}_1& U_2\end{array}\right],U_1\in {\mathbb{R}}^{m\times n},U_2\in {\mathbb{R}}^{m\times (m-n)},\end{array}$$

and it follows from the orthogonal property of U that

$$\begin{array}{c}\hfill \begin{array}{ccc}{U}_1^T{U}_1=I_n,\hfill & & U_2^T{U}_2=I_{m-n},\hfill \\ U_1^T{U}_2=0_{n,m-n},\hfill & & U_2^T{U}_1=0_{m-n,n},\hfill \end{array}\end{array}$$

(13)

and

$$\begin{array}{c}\hfill U_1{U}_1^T+U_2{U}_2^T=I_m.\end{array}$$

(14)

It follows from (12) that c can be written as

$$\begin{array}{c}\hfill c=\left[\begin{array}{c}{\overline{c}}_1\\ {\overline{c}}_2\end{array}\right]=\left[\begin{array}{c}{U}_1^T\\ U_2^T\end{array}\right]b,{\overline{c}}_1\in {\mathbb{R}}^n,{\overline{c}}_2\in {\mathbb{R}}^{m-n},\end{array}$$

where ${\overline{c}}_1$ and ${\overline{c}}_2$ are the components of b that lie in, respectively, the spaces spanned by the columns of $U_1$ and $U_2$. It follows from (11) that

$$\begin{array}{c}\hfill e^2=\frac{{\sum }_{i=n+1}^m{c}_i^2}{{\sum }_{i=1}^m{c}_i^2}=\frac{{∥{\overline{c}}_2∥}^2}{{∥c∥}^2}=\frac{{∥{\overline{c}}_2∥}^2}{{∥{\overline{c}}_1∥}^2+{∥{\overline{c}}_2∥}^2}=1-{cos}^2\theta ,\end{array}$$

(15)

where

$$\begin{array}{c}\hfill cos\theta =\frac{∥{\overline{c}}_1∥}{∥c∥},\end{array}$$

(16)

and thus the error in $x_0$ is a measure of the proportion of b that lies in the space that is orthogonal to $\mathcal{R}\left(A\right)$, and a large value of $\theta$ leads to a large error in $x_0$. The significance of $\theta$ is considered in Section 4 and it is shown that it allows a geometric interpretation of the conditions that define the minimum and maximum values of e.

4. The Angle between $\mathit{b}$ and $\mathcal{R}\left(\mathit{A}\right)$

It is shown in this section that $\theta$, where $cos\theta$ is defined in (16), is the angle between b and its component that lies in the column space $\mathcal{R}\left(A\right)$ of A. This interpretation follows from the decomposition of b into a component $r_1\in \mathcal{R}\left(A\right)$ and a component $r_2\in \mathcal{R}{\left(A\right)}^{\perp }$,

$$\begin{array}{c}\hfill b=r_1+r_2.\end{array}$$

Since $r_1\in \mathcal{R}\left(A\right)$ and $r_2\in \mathcal{R}{\left(A\right)}^{\perp }$, there exist vectors $t_1$ and $t_2$ such that

$$\begin{array}{c}\hfill r_1=U_1{t}_1\mathrm{and}{r}_2=U_2{t}_2,\end{array}$$

because the columns of $U_1$ and $U_2$ form orthonormal bases for $\mathcal{R}\left(A\right)$ and $\mathcal{R}{\left(A\right)}^{\perp }$, respectively, and thus

$$\begin{array}{c}\hfill b=U_1{t}_1+U_2{t}_2.\end{array}$$

It follows from (13) that premultiplication of this equation by $U_1^T$ and then $U_2^T$ yields

$$\begin{array}{c}\hfill b=U_1{U}_1^{T}b+U_2{U}_2^{T}b,\end{array}$$

where the first and second terms on the right are, respectively, the orthogonal projections of b onto $\mathcal{R}\left(A\right)$ and $\mathcal{R}{\left(A\right)}^{\perp }$. It follows from (13) that the second term is orthogonal to the first term, which confirms that $r_1$ is orthogonal to $r_2$. Also, it follows from (14) that

$$\begin{array}{c}\hfill {∥U_1{U}_1^{T}b∥}^2+{∥U_2{U}_2^{T}b∥}^2={∥b∥}^2,\end{array}$$

where

$$\begin{array}{c}\hfill ∥U_1{U}_1^{T}b∥=∥U_1^{T}b∥=∥{\overline{c}}_1∥,\end{array}$$

and thus $cos\theta$ is equal to the ratio of the magnitude of the orthogonal projection of b onto $\mathcal{R}\left(A\right)$ to the magnitude of b. Equation (15) shows that there is a simple relationship between the error in $x_0$, and the angle between b and its component that lies in $\mathcal{R}\left(A\right)$. The limiting cases are:

• $b\in \mathcal{R}\left(A\right)$: It follows that $cos\theta =1$ and $U_2=0$, and thus $e=0$.
• $b\in \mathcal{R}{\left(A\right)}^{\perp }$: It follows that $cos\theta =0$ and $U_1=0$, and thus $e=1$.

The ratio $cos\theta$ provides a geometric interpretation of the error in the solution of the $\mathrm{LS}$ problem, and it is shown in Section 5 that $cos\theta$ also arises in the expressions for the condition numbers of $x_0$.

5. The Numerical Condition of the Optical Flow

This section considers the numerical condition of the optical flow $x_0$, the solution of the $\mathrm{LS}$ problem. The condition number of A, $\kappa \left(A\right)=\raisebox{1ex}{${\sigma }_1$}\!\left/ \!\raisebox{-1ex}{${\sigma }_2$}\right.$, where ${\sigma }_1>{\sigma }_2$, is frequently used as a measure of the stability of $x_0$, but it requires qualification:

• The condition number $\kappa \left(A\right)$ is a function of A and it is independent of b, but $x_0=A^{†}b$ is a function of A and b. This independence of b follows because $\kappa \left(A\right)$ is the upper bound of the ratio of the relative error in $x_0=A^{†}b$ to the relative error in b with respect to all vectors $b\in \mathcal{R}\left(A\right)$ and perturbations $\delta b\in {\mathbb{R}}^m$. The vector b is specified in a given problem and thus a better measure of the stability of $x_0$ is an upper bound of the ratio $\Delta$ of the relative error in $x_0$ to the relative error in b with respect to all perturbations $\delta b$, for the given vector b. The condition number $\kappa \left(A\right)$ may therefore overestimate the true value of $\Delta$ by many orders of magnitude if $b\in \mathcal{R}\left(A\right)$.
• The property $b\in \mathcal{R}\left(A\right)$ is not, in general, satisfied in the $\mathrm{LK}$ method and thus a measure of the stability of $x_0$ with respect to perturbations in b that is appropriate when $b\notin \mathcal{R}\left(A\right)$ must be developed.

Although these points show that $\kappa \left(A\right)$ as a measure of the stability of $x_0$ has disadvantages, it has been used to assess the computational reliability of $x_0$ ([§II.A] [14]), ([§4.3.2] [16]), ([§2] [10]). These disadvantages of $\kappa \left(A\right)$ are overcome by the effective condition number of $x_0$, which is a function of A and b, and therefore a more accurate measure of the stability of $x_0$. This condition number is considered in Section 5.1.

The effective condition number is a linear condition number because it is assumed that the entries of A (the spatial derivatives of the intensity) are not subject to error, and only errors in b (the temporal derivatives of the intensity) are subject to error. A more realistic condition number of the optical flow requires that errors in the spatial and temporal derivatives of the intensity be considered. This inclusion leads to a non-linear extension of the effective condition number, and it is considered in Section 5.2.

5.1. The Effective Condition Number

This section considers the effective condition number of the solution $x_0$ of the $\mathrm{LS}$ problem, which must be solved many times in the $\mathrm{LK}$ method. It is assumed for generality that the coefficient matrix A is of order $m\times n$ where $m\ge n$, rather than $m\times 2$.

The relative errors of $x_0$ and b are, respectively,

$$\begin{array}{c}\hfill \Delta x_0=\frac{∥\delta x_0∥}{∥x_0∥}and\Delta b=\frac{∥\delta b∥}{∥b∥},\end{array}$$

(17)

and the effective condition number of $x_0$ is defined as the maximum value of the ratio of $\Delta x_0$ to $\Delta b$ for the given vector b, with respect to all perturbations $\delta b$ in b,

$$\begin{array}{c}\hfill \eta (A,b)=\underset{\delta b\in {\mathbb{R}}^m}{max}\frac{\Delta x_0}{\Delta b}.\end{array}$$

Theorem 1.

The effective condition number $\eta (A,b)$ of $x_0$ is

$$\begin{array}{c}\hfill \eta (A,b)=\frac{∥A^{†}∥∥b∥}{∥x_0∥}=\frac{∥A^{†}∥∥b∥}{∥A^{†}b∥},\end{array}$$

(18)

and it follows from (12) that

$$\begin{array}{c}\hfill \eta (A,b)=\frac{∥c∥}{{\sigma }_n∥{\Sigma }^{†}c∥}.\end{array}$$

(19)

Proof. 

It follows from (7) and (17) that

$$\begin{array}{c}\hfill ∥\delta x_0∥=∥A^{†}\delta b∥\le ∥A^{†}∥∥\delta b∥=∥A^{†}∥∥b∥\Delta b.\end{array}$$

(20)

The division of this inequality by $∥x_0∥=∥A^{†}b∥$ yields (18), from which (19) follows. □

The condition that must be satisfied for $\kappa \left(A\right)$ to be a measure of the stability of $x_0$ follows from (20) because $∥b-A{x}_0∥=0$ if and only if $b\in \mathcal{R}\left(A\right)$. This condition is necessarily satisfied by all vectors b if A is square, but by some vectors b, not all vectors b, if $m>n$. If $b\in \mathcal{R}\left(A\right)$, then $b=A{x}_0$ and thus $∥b∥\le ∥A∥∥x_0∥$, and (20) yields

$$\begin{array}{c}\hfill ∥\delta x_0∥\le ∥A^{†}∥∥b∥\Delta b\le ∥A^{†}∥∥A∥∥x_0∥\Delta b.\end{array}$$

This leads to the definition of the condition number $\kappa \left(A\right)$ of $A\in {\mathbb{R}}^{m\times n}$, $m\ge n$, as the upper bound of the ratio of the relative error of $x_0$ to the relative error of b if and only if $b\in \mathcal{R}\left(A\right)$,

$$\begin{array}{c}\hfill \kappa \left(A\right)=\underset{\delta b\in {\mathbb{R}}^m,b\in \mathcal{R}\left(A\right)}{max}\frac{\Delta x_0}{\Delta b}=∥A^{†}∥∥A∥=\frac{{\sigma }_1}{{\sigma }_n}.\end{array}$$

Furthermore, it follows from (12) and (19) that

$$\begin{array}{c}\hfill \underset{b\in \mathcal{R}\left(A\right)}{max}\eta (A,b)=\underset{c\in \mathcal{R}\left(A\right)}{max}\frac{∥c∥}{{\sigma }_n∥{\Sigma }^{†}c∥}=\frac{{\sigma }_1}{{\sigma }_n},\end{array}$$

and thus $\kappa \left(A\right)$ is the maximum value of $\eta (A,b)$ with respect to all vectors $b\in \mathcal{R}\left(A\right)$.

The superiority of the effective condition number $\eta (A,b)$ with respect to the condition number $\kappa \left(A\right)$ as a measure of the stability of $x_0$ follows because $\eta (A,b)$ is a function of A and b, but $\kappa \left(A\right)$ is a function of A only. The example of regression in ([§3] [23]) shows the difference between these condition numbers. The dependence of $\eta (A,b)$ on A and b is clearly advantageous, but $\eta (A,b)$ must be used with care because it follows from (18) that the denominator contains the term $∥x_0∥$. It therefore follows that $\eta (A,b)$ is ill conditioned if $x_0$ is ill conditioned, and it is well conditioned if $x_0$ is well conditioned. It is shown in ([§4] [24]) that $\eta (A,b)$, and therefore $x_0$, are ill conditioned if the discrete Picard condition is satisfied ([§4.5] [25]). This condition states that if the constants $\left|c_i\right|$ decay to zero faster than the singular values ${\sigma }_i$ decay to zero, that is,

$$\begin{array}{c}\hfill \frac{\left|c_i\right|}{{\sigma }_i}\to 0\mathrm{as}i\to n,\end{array}$$

(21)

then $x_0$ is ill conditioned. This condition can be derived from $x_0$,

$$\begin{array}{c}\hfill x_0=A^{†}b=\sum _{i=1}^n\left(\frac{c_i}{{\sigma }_i}\right)v_i,\end{array}$$

(22)

where $v_i$ is the ith column of V. If b is perturbed to $b+\delta b$, then $x_0$ is perturbed to

$$\begin{array}{c}\hfill x_0+\delta x_0=\sum _{i=1}^n\left(\frac{c_i+\delta c_i}{{\sigma }_i}\right)v_i,\end{array}$$

(23)

and there exists a constant s such that the perturbations $\left|\delta c_i\right|\approx ϵ,i=1,\cdots ,n$, satisfy

$$\begin{array}{c}\hfill \begin{array}{ccccc}\left|c_i\right|\hfill & >& ϵ,& & i=1,\cdots ,s,\hfill \\ \left|c_i\right|\hfill & <& ϵ,& & i=s+1,\cdots ,n.\hfill \end{array}\end{array}$$

These inequalities follow because the perturbations $\left|\delta c_i\right|$ are approximately constant but the coefficients $\left|c_i\right|$ decay to zero because the discrete Picard condition (21) is satisfied, and thus from (23),

$$\begin{array}{c}\hfill x_0+\delta x_0=\sum _{i=1}^n\left(\frac{c_i+\delta c_i}{{\sigma }_i}\right)v_i\approx \sum _{i=1}^s\left(\frac{c_i}{{\sigma }_i}\right)v_i+\sum _{i=s+1}^n\left(\frac{\delta c_i}{{\sigma }_i}\right)v_i.\end{array}$$

Since the discrete Picard condition is satisfied, $x_0$ can be approximated with a small error by only considering the first s singular values in (22),

$$\begin{array}{c}\hfill x_0=A^{†}b\approx \sum _{i=1}^s\left(\frac{c_i}{{\sigma }_i}\right)v_i,\end{array}$$

and thus

$$\begin{array}{c}\hfill \delta x_0\approx \sum _{i=s+1}^n\left(\frac{\delta c_i}{{\sigma }_i}\right)v_i.\end{array}$$

It follows that

$$\begin{array}{c}\hfill ∥\delta x_0∥\approx ϵ{\left(\sum _{i=s+1}^n\frac{1}{{\sigma }_i^2}\right)}^{\frac{1}{2}}\approx \frac{ϵ}{{\sigma }_n},\end{array}$$

and thus the norm of the perturbation $\delta x_0$ is approximately proportional to the magnitude of the noise and inversely proportional to the smallest singular value of A, from which it follows that $x_0$ is ill conditioned if the discrete Picard condition (21) is satisfied. It is seen that $x_0$ is ill conditioned if it is dominated by the large singular values of A, in which case $\delta x_0$ is dominated by the small singular values of A.

The satisfaction of the discrete Picard condition is necessary for the application of Tikhonov regularisation because it guarantees that the regularisation error is small and the regularised form of $x_0$ is numerically stable ([§5] [23]). Tikhonov regularisation is used for the removal of blur from images because experiments have shown that many images satisfy the discrete Picard condition [26]. This condition forms the prior information that guarantees that Tikhonov regularisation is effective for image deblurring, and more generally, the ill conditioned property of $\eta (A,b)$ limits its practical use if prior information is not available. The condition number $\eta (A,b)$ is, however, an important condition measure because, as shown above, it defines the conditions between A and b for which $x_0$ is well conditioned, and the conditions for which $x_0$ is ill conditioned.

It has been shown that $\eta (A,b)$ is ill conditioned if $x_0$ is ill conditioned, but its upper bound is numerically stable, and an expression for this bound is developed in Theorem 2. This bound includes the term $cos\theta$ that is defined in (16), and thus the error in the solution $x_0$ of the $\mathrm{LS}$ problem is related to the numerical condition of $x_0$.

Theorem 2.

If $A\in {\mathbb{R}}^{m\times n}$, $m\ge n$, then the effective condition number (19) of the solution $x_0$ of the $\mathit{LS}$ problem is

$$\begin{array}{c}\hfill \eta (A,b)=\frac{1}{{\sigma }_n}{\left(\frac{{\sum }_{i=1}^m{c}_i^2}{{\sum }_{i=1}^n{\left(\frac{c_i}{{\sigma }_i}\right)}^2}\right)}^{\frac{1}{2}}=\left(\frac{{\sigma }_1}{{\sigma }_n}\right){\left(\frac{{\sum }_{i=1}^m{c}_i^2}{{\sum }_{i=1}^n{\left(\frac{{\sigma }_1}{{\sigma }_i}\right)}^2{c}_i^2}\right)}^{\frac{1}{2}},\end{array}$$

where $c={\left\{c_i\right\}}_{i=1}^m$ is defined in (12), and thus $\eta (A,b)$ satisfies

$$\begin{array}{c}\hfill \eta (A,b)\le \left\{\begin{array}{ccccccc}\kappa \left(A\right)\hfill & & & if& & & m=n,\hfill \\ \frac{\kappa \left(A\right)}{cos\theta }\hfill & & & if& & & m>n.\hfill \end{array}\right.\end{array}$$

(24)

Proof. 

If $m=n$, (19) yields

$$\begin{array}{c}\hfill \eta (A,b)=\left(\frac{{\sigma }_1}{{\sigma }_n}\right){\left(\frac{{\sum }_{i=1}^n{c}_i^2}{{\sum }_{i=1}^n{\left(\frac{{\sigma }_1}{{\sigma }_i}\right)}^2{c}_i^2}\right)}^{\frac{1}{2}}\le \frac{{\sigma }_1}{{\sigma }_n}=\kappa \left(A\right),\end{array}$$

and the result (24) for $m=n$ follows.

Consider now the situation $m>n$,

$$\begin{array}{c}\hfill \eta (A,b)=\left(\frac{{\sigma }_1}{{\sigma }_n}\right){\left(\frac{{\sum }_{i=1}^m{c}_i^2}{{\sum }_{i=1}^n{\left(\frac{{\sigma }_1}{{\sigma }_i}\right)}^2{c}_i^2}\right)}^{\frac{1}{2}}\le \left(\frac{{\sigma }_1}{{\sigma }_n}\right){\left(\frac{{\sum }_{i=1}^m{c}_i^2}{{\sum }_{i=1}^n{c}_i^2}\right)}^{\frac{1}{2}}=\kappa \left(A\right)\frac{∥c∥}{∥{\overline{c}}_1∥},\end{array}$$

and the result (24) for $m>n$ follows from (16). □

Equation (24) shows that the condition number of $x_0$ increases rapidly as $cos\theta \to 0$. For example, if a window of order $7\times 7$ is used, then $A\in {\mathbb{R}}^{49\times 2}$ and the error in $x_0$ is equal to zero if and only if b lies in the two dimensional subspace ${\mathbb{R}}_{colA}^2$ spanned by the columns of A, in the space ${\mathbb{R}}^{49}$, and the error in $x_0$ is due to the component of b that lies in the space that is orthogonal to ${\mathbb{R}}_{colA}^2$. The condition number of $x_0$ is large if $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$cos\theta $}\right.\gg 1$, even if $\raisebox{1ex}{${\sigma }_1$}\!\left/ \!\raisebox{-1ex}{${\sigma }_2$}\right.=O\left(1\right)$.

Example 1.

Consider the matrix A and vector b,

$$\begin{array}{c}\hfill A=\left[\begin{array}{cc}\hfill 1& \hfill 0\\ \hfill 0& \hfill 1\\ \hfill 1& \hfill -1\end{array}\right]andb=\left[\begin{array}{c}\hfill -1\\ \hfill 1\\ \hfill 1\end{array}\right].\end{array}$$

The singular values of A are ${\sigma }_1=1$ and ${\sigma }_2=\sqrt{3}$, and thus $\kappa \left(A\right)=\sqrt{3}$. It may therefore be thought that the solution $x_0$ of the $\mathit{LS}$ problem is well conditioned but this is incorrect because b is orthogonal to $\mathcal{R}\left(A\right)$,

$$\begin{array}{c}\hfill b^{T}A=\left[\begin{array}{ccc}0& 0& 0\end{array}\right].\end{array}$$

It follows that $cos\theta =0$ and $x_0={\left[\begin{array}{cc}0& 0\end{array}\right]}^T$.

The left singular matrix U of $A=U\Sigma V^T$, and the matrices $U_1$ and $U_2$ are

$$\begin{array}{c}\hfill U=\left[\begin{array}{ccccc}\hfill \frac{1}{\sqrt{6}}& & \hfill \frac{1}{\sqrt{2}}& & \hfill \frac{-1}{\sqrt{3}}\\ \hfill \frac{-1}{\sqrt{6}}& & \hfill \frac{1}{\sqrt{2}}& & \hfill \frac{1}{\sqrt{3}}\\ \hfill \frac{2}{\sqrt{6}}& & \hfill 0& & \hfill \frac{1}{\sqrt{3}}\end{array}\right],U_1=\left[\begin{array}{ccc}\hfill \frac{1}{\sqrt{6}}& & \hfill \frac{1}{\sqrt{2}}\\ \hfill \frac{-1}{\sqrt{6}}& & \hfill \frac{1}{\sqrt{2}}\\ \hfill \frac{2}{\sqrt{6}}& & \hfill 0\end{array}\right]and{U}_2=\left[\begin{array}{c}\hfill \frac{-1}{\sqrt{3}}\\ \hfill \frac{1}{\sqrt{3}}\\ \hfill \frac{1}{\sqrt{3}}\end{array}\right],\end{array}$$

and thus

$$\begin{array}{c}\hfill U_1^{T}b=\left[\begin{array}{c}\hfill 0\\ \hfill 0\end{array}\right]and{U}_1{U}_1^{T}b=\left[\begin{array}{c}\hfill 0\\ \hfill 0\\ \hfill 0\end{array}\right],\end{array}$$

which confirm that b is orthogonal to $\mathcal{R}\left(A\right)$.

5.2. A Non-Linear Condition Number

This section extends the linear effective condition number considered in Section 5.1 by including perturbations in the coefficient matrix A, such that the perturbed solution of the $\mathrm{LS}$ problem is

$$\begin{array}{c}\hfill x_0+\delta x_0={\left(A+\delta A\right)}^{†}(b+\delta b).\end{array}$$

(25)

An expression for the condition number of $x_0$ when A and b are perturbed is developed in Theorem 3.

Theorem 3.

The condition number of $x_0$ when A and b are perturbed by, respectively, $\delta A$ and $\delta b$ such that

$$\begin{array}{c}\hfill ∥\delta b∥\le ϵ∥b∥and∥\delta A∥\le ϵ∥A∥,\end{array}$$

(26)

is

$$\begin{array}{ccc}\hfill \rho (A,b)& =& \underset{∥\delta b∥\le ϵ∥b∥,∥\delta A∥\le ϵ∥A∥}{max}\frac{1}{ϵ}\frac{∥\delta x_0∥}{∥x_0∥}\hfill \\ & =& \frac{∥A^{†}∥∥b∥+∥A∥∥{\left(A^{T}A\right)}^{-1}∥∥\left(I-A{A}^{†}\right)b∥}{∥x_0∥}+∥A∥∥A^{†}∥\hfill \\ & =& \eta (A,b)+\kappa \left(A\right)+\frac{\kappa \left(A\right)∥\left(I-A{A}^{†}\right)b∥}{{\sigma }_n∥x_0∥}.\hfill \end{array}$$

(27)

Proof. 

The pseudo-inverse of $A+\delta A$ is

$$\begin{array}{c}\hfill {(A+\delta A)}^{†}={\left({\left(A+\delta A\right)}^T\left(A+\delta A\right)\right)}^{-1}{\left(A+\delta A\right)}^T,\end{array}$$

where, to first order,

$$\begin{array}{ccc}\hfill {\left({\left(A+\delta A\right)}^T\left(A+\delta A\right)\right)}^{-1}& =& {\left(I+A^{†}\delta A+{\left(A^{T}A\right)}^{-1}\delta A^{T}A\right)}^{-1}{\left(A^{T}A\right)}^{-1}\hfill \\ & =& \left(I-A^{†}\delta A-{\left(A^{T}A\right)}^{-1}\delta A^{T}A\right){\left(A^{T}A\right)}^{-1},\hfill \end{array}$$

and thus

$$\begin{array}{ccc}\hfill {(A+\delta A)}^{†}& =& {\left({\left(A+\delta A\right)}^T\left(A+\delta A\right)\right)}^{-1}{(A+\delta A)}^T\hfill \\ & =& \left(I-A^{†}\delta A-{\left(A^{T}A\right)}^{-1}\delta A^{T}A\right){\left(A^{T}A\right)}^{-1}{(A+\delta A)}^T\hfill \\ & =& A^{†}+{\left(A^{T}A\right)}^{-1}\delta A^T-A^{†}\delta A{A}^{†}-{\left(A^{T}A\right)}^{-1}\delta A^{T}A{A}^{†},\hfill \end{array}$$

to first order. It therefore follows from (25) that

$$\begin{array}{c}\hfill x_0+\delta x_0=\left(A^{†}+{\left(A^{T}A\right)}^{-1}\delta A^T-A^{†}\delta A{A}^{†}-{\left(A^{T}A\right)}^{-1}\delta A^{T}A{A}^{†}\right)(b+\delta b),\end{array}$$

and thus

$$\begin{array}{c}\hfill \delta x_0=A^{†}\delta b-A^{†}\delta A{x}_0+{\left(A^{T}A\right)}^{-1}\delta A^T\left(I-A{A}^{†}\right)b.\end{array}$$

Hence,

$$\begin{array}{ccc}\hfill ∥\delta x_0∥& \le & ∥A^{†}∥∥\delta b∥+∥A^{†}∥∥\delta A∥∥x_0∥+∥{\left(A^{T}A\right)}^{-1}∥∥\delta A∥∥\left(I-A{A}^{†}\right)b∥,\hfill \end{array}$$

and the substitution of (18) and (26) into this expression yields (27). □

An upper bound for the last term in (27) allows $\rho (A,b)$ to be expressed in terms of $\eta (A,b)$ and $\kappa \left(A\right)$. In particular,

$$\begin{array}{c}\hfill ∥\left(I-A{A}^{†}\right)b∥=∥\left[\begin{array}{cc}0& 0\\ 0& I_{m-n}\end{array}\right]U^{T}b∥\le ∥b∥,\end{array}$$

and thus (27) becomes

$$\begin{array}{ccc}\hfill \rho (A,b)& =& \eta (A,b)+\kappa \left(A\right)+\frac{\kappa \left(A\right)∥\left(I-A{A}^{†}\right)b∥}{{\sigma }_n∥x_0∥}\hfill \\ & \le & \eta (A,b)+\kappa \left(A\right)+\frac{\kappa \left(A\right)∥b∥}{{\sigma }_n∥A^{†}b∥}\hfill \\ & =& \eta (A,b)+\kappa \left(A\right)+\kappa \left(A\right)\eta (A,b),\hfill \end{array}$$

(28)

and thus the inclusion of errors in the spatial derivatives of the intensity is significant because of the term $\kappa \left(A\right)\eta (A,b)$.

It was stated in Section 5.1 that the effective condition number $\eta (A,b)$ is ill conditioned if the discrete Picard condition is satisfied, in which case (27) and its upper bound (28) cannot be computed reliably. This problem is addressed by using the upper bound (24) for the value of $\eta (A,b)$ in (28). The condition number (27) is, however, useful because it shows that the inclusion of perturbations in the spatial derivatives of the intensity causes a significant increase in the condition number of the optical flow.

6. Examples

This section contains two examples that show the errors and condition numbers of the optical flow at each level of the pyramid. The Middlebury data set [1] is used in Examples 2 and 3, and the following numerical measures are considered:

• The condition number $\kappa \left(A\right)$, the upper bound of the effective condition number, $\raisebox{1ex}{$\kappa \left(A\right)$}\!\left/ \!\raisebox{-1ex}{$cos\theta $}\right.$, and the non-linear condition number $\rho (A,b)$.
• The angle $\theta$ between b and its component that lies in $\mathcal{R}\left(A\right)$.
• The error e in the solution $x_0$ of the $\mathrm{LS}$ problem.

The examples use the images shown in Figure 2 [1].

Example 2.

The five measures listed above (the three condition numbers, the angle θ and the error e in $x_0$) were computed in each window for windows of side lengths $2s+1$, $s=2,3,\dots ,8,9,$ at each level of a pyramid that has four levels, for the images in Figure 2. The mean and standard deviation of ${log}_{10}\kappa \left(A\right)$, the upper bound ${log}_{10}\raisebox{1ex}{$\kappa \left(A\right)$}\!\left/ \!\raisebox{-1ex}{$cos\theta $}\right.$ of the effective condition number, and the non-linear condition number ${log}_{10}\rho (A,b)$, at each level of the pyramid and for each window size, are shown in Figure 3. The figures in the left column show that, as expected, the mean of $\rho (A,b)$ at each of the 32 points in the grid is larger than the corresponding values of $\kappa \left(A\right)$ and $\raisebox{1ex}{$\kappa \left(A\right)$}\!\left/ \!\raisebox{-1ex}{$cos\theta $}\right.$. Also, the means of $\rho (A,b)$ and $\raisebox{1ex}{$\kappa \left(A\right)$}\!\left/ \!\raisebox{-1ex}{$cos\theta $}\right.$ increase, but not monotonically, as the pyramid is descended from the coarsest image (level 4) to the finest image (level 1), and they decrease as the window size increases, at each level of the pyramid. The variation of the mean of each condition number, Figure 3 (left), is similar to the variation of its standard deviation, Figure 3 (right).

The mean and standard deviation of the angle θ, where $cos\theta$ is defined in (16), and the error ${log}_{10}e$, where e is defined in (11), are shown in Figure 4. It is seen that the mean of θ is approximately constant for windows of side length 7 pixels or more at levels 2, 3 and 4 of the pyramid. Figure 4 (left) shows that the variation of the means of θ and ${log}_{10}e$ are very similar, but Figure 4 (right) show that, for each measure, the standard deviation is significant, which implies that there is considerable variation of each measure.

Example 3.

The $\mathit{LK}$ method with a pyramid of three levels was applied to the images in Figure 2, with a window of size $15\times 15$. The condition number ${log}_{10}\kappa \left(A\right)$, upper bound ${log}_{10}\raisebox{1ex}{$\kappa \left(A\right)$}\!\left/ \!\raisebox{-1ex}{$cos\theta $}\right.$ of the effective condition number, the non-linear condition number ${log}_{10}\rho (A,b)$, angle θ and error ${log}_{10}e$ were computed for every $\mathit{LS}$ problem in the windows for each level of the pyramid. Figure 5 shows the axes and size of the image at each level of the pyramid for Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10. These figures show that the quantities above increase, and the standard deviation of each of these quantities increases, as the pyramid is descended. The figures are therefore consistent with Figure 3 and Figure 4, and it follows that the computational reliability of the optical flow decreases as the pyramid is descended.

Figure 9 shows that the values of the local maxima of the angle θ increase significantly as the pyramid is descended. In particular there are many $\mathrm{LS}$ problems for which θ is about 80 degrees, from which it follows that the errors in the solutions of these problems are large. This result is consistent with the plots of the condition numbers in Figure 7 and Figure 8. Furthermore, it follows from (15) that large values of θ are associated with large errors e, which is confirmed in Figure 10.

7. Summary and Conclusions

The work in this paper considers the numerical properties of the solution of the $\mathrm{LS}$ problem because of their consequences on the computed optical flow. These consequences are compounded by the large number of $\mathrm{LS}$ problems that must be solved in the $\mathrm{LK}$ method, where the solutions of the $\mathrm{LS}$ problems at level l of the pyramid form the data, after upsampling, for the $\mathrm{LS}$ problems at level $l-1$ of the pyramid. It has been shown that the large values of the condition numbers $\eta (A,b)$ and $\rho (A,b)$, and error e, of the solution $x_0$ of the $\mathrm{LS}$ problem arise because a significant component of b lies in the space that is orthogonal to $\mathcal{R}\left(A\right)$, and thus the terms $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$cos\theta $}\right.$ and $cos\theta$ in, respectively, the expressions for the condition numbers and error are significant. It has been shown that the error and condition numbers of the $\mathrm{LS}$ problem increase as the pyramid is descended from the top level (coarsest image) to the bottom level (finest image) and thus the fidelity of the computed optical flow decreases as the pyramid is descended. There have been several studies that consider the errors associated with the $\mathrm{LK}$ method, but the new aspect of the work described in this paper is the extensive use of computational linear algebra, specifically, condition estimation and error analysis of the solution of the $\mathrm{LS}$ problem.

The large values of the condition numbers and errors follow from the assumption that the optical flow in each window at each level of the pyramid is constant, that is, the motion in each window can be represented by a translation. It follows that an improved form of the $\mathrm{LK}$ method requires that a richer class of motions be allowed in each window, and this is achieved by the introduction of warping functions in the $\mathrm{LK}$ method [17]. These warping functions are defined by parametric models, for example, affine and projective models, and the values of the parameters of these models are computed in the optimisation procedure in an extended form of the $\mathrm{LK}$ method.