Dense Parameter Fields from Total Least Squares - CiteSeerX
Dense Parameter Fields from Total Least Squares Hagen Spies1 and Christoph S. Garbe2 1
ICG-III: Phytosphere, Research Center J¨ ulich, 52425 J¨ ulich, Germany, [email protected] 2 Interdisciplinary Center for Scientific Computing, University of Heidelberg, INF 368, 69120 Heidelberg, Germany, [email protected]
Abstract. A method for the interpolation of parameter fields estimated by total least squares is presented. This is applied to the study of dynamic processes where the motion and further values such as divergence or brightness changes are parameterised in a partial differential equation. For the regularisation we introduce a constraint that restricts the solution only in the subspace determined by the total least squares procedure. The performance is illustrated on both synthetic and real test data.
1
Introduction
The parameters governing the changes in an image sequence can be estimated locally by means of a total least squares technique [1; 2]. For the simple case of movement alone, i.e. optical flow, it is known that in areas with one dominant gradient direction only the component of the velocity in that direction can be computed. This is the well known aperture problem [3]. When there are more parameters to be estimated there might be even more such linear dependencies between them. We show how such cases can be detected and how an appropriate minimum norm solution can be obtained. In many applications we are nevertheless interested in dense parameter fields and thus there is the need for an interpolation procedure. Clearly it is desirable to utilise the available information from the minimum norm solutions. Towards this end we present a regularisation framework that restricts the interpolated solution only in the subspace in which a solution is available. This is in contrast to variational methods that directly seek to minimise the constraint over the entire data set [4; 5]. The main idea of the presented regularisation has been introduced for the special case of 3D velocity (range flow) estimation in [6]. Here this is extended and applied to more general dynamic processes.
2
Dynamic Processes in Image Sequences
In scientific image processing one often encounters a linear equation for a set of model parameters m at each data point (pixel) of type: aT m = b; a, m ∈ IRm .
An ordinary least squares estimate implies that only b is erroneous while a is exact. This assumption certainly does not hold for the cases discussed below where a contains derivatives. Therefore we introduce a new parameter vector p of dimension n = m + 1 and rewrite the constraint as: dT p = 0 with d = (aT , −b)T
and p = (mT , 1)T .
(1)
This concept is very general in the sense that the parameters of any dynamic process that can be modeled by a linear partial differential equation falls into this scheme [2; 7]. Prominent physical examples are source terms, relaxation and diffusion processes. If there is no such model available a Taylor series expansion up to the desired order may be used [7]. A few examples that are used in the later experiments (Sect. 5) are: Optical flow. The assumption that all temporal variation in the image brightdy dx ness g(x, y, t) is caused by movements implies [4]: dg dt = gx dt + gy dt + gt = 0. With indices denoting partial derivatives. This brightness change constraint equation forms the basis of many algorithms to compute the optical flow dy T [u, v]T = [ dx dt , dt ] [3; 1]. This is readily expressed in the form of (1) using T d = [gx gy gt ] and p = [u v 1]T . Source terms. Source terms cause first-order temporal changes in the image seT quence proportional to a source strength q: ( dg dt = q) yielding d = [gx gy 1 gt ] T and p = [u v q 1] . Source terms occur for example in infrared images if an object is heated or cooled or if the global illumination of a scene changes [8]. The same equation is obtained in the study of the three-dimensional movement of surfaces from sequences of depth maps: zx u + zy v + zt = w. Here [u v w]T denotes the 3D-movement also called range flow [6]. Divergent velocity. The proposed concept can also model spatial variations in the velocity field, e.g. an affine velocity field [1]. If we are interested in a particular aspect of the flow field, such as its divergence or rotation, we can only include this in our estimation. For divergence d this reads as: d = [gx gy (xgx +ygy ) gt ]T ; p = [u v d 1]T . Here x, y denote local coordinates with respect to the computed location. Examples include camera movement along the optical axis or local growth rates in biological experiments.
3
TLS Estimation
The presented estimation framework extends previous work on total least squares optical flow estimation [2]. For an estimation we pool the constraints (1) over a local neighbourhood and assume the parameters to be constant within: Z ∞ ˆ = arg min ˆ )2 dx0 dt0 subject to p ˆT p ˆ = 1 . (2) p w(x − x0 , t − t0 )(dT p −∞
Here w(x − x0 , t − t0 ) is a weighting function that defines the spatio-temporal neighbourhood where the parameters are estimated. In order to avoid the trivial
a
b
c
d
Fig. 1: Moving square example: a image data, b type of flow (black: trace to small, light grey: linear dependency and dark grey: full estimate possible), c normal flow and d full flow.
ˆ to be normalised. Minimisation leads solution we require the solution vector p to an eigenvalue equation of the extended structure tensor J : Z ∞ ˆ = λˆ Jp p with J (x) = w (x − x0 , t − t0 ) (d dT ) dx0 dt0 . (3) −∞
It follows that the solution to (2) is given by the nullspace of J [9]. Section 3.2 describes how to recover the parameters. For an unbiased and optimal TLS estimate the errors in all data terms need to be independent random variables with zero mean and equal variance [9; 10]. This implies a scaling with respect to the absolute error values. 3.1
Type and Confidence Values
The spectrum of J can be used to detect linear dependencies in the data. In the following we assume the eigenvalues of J to be sorted in descending order: λ1 ≥ λ2 ≥ . . . ≥ λn−1 ≥ λn . The trace of J is used to determine if there is enough variation (texture) in the data for the computation to make sense. Thus we only proceed where tr(J ) > τ1 . The smallest eigenvalue directly yields the residual of the parameter fit. Thus we reject unreliable estimates if λn > τ2 , with a threshold τ2 corresponding to the noise level. This occurs at discontinuities where a single set of parameters is not sufficient or when no coherent motion is present. A small λn does not ensure that all the parameters can be estimated independently from each other. More than one eigenvalue may be close to zero (< τ2 ) and we can no longer uniquely pick a solution. Any vector in the nullspace of J is a possible solution, see Sect. 3.2. A type measure that indicates how well the dimensionality of the nullspace of J has been determined can be obtained by examining how much the first relevant eigenvalue λq is above the threshold: 2 λ q − τ2 ωt = , λ1 ≥ . . . ≥ λq > λq+1 ≈ . . . ≈ λn < τ2 . (4) λq For a simple optical flow example of a moving square (Fig. 1) a map indicating which type of flow can be computed is given in Fig. 1b.
3.2
Minimum Norm Solutions
As stated above any vector in the nullspace of J is a possible solution. A sensible choice is that resulting in an estimate with minimal norm. This solution is given ˆi (i > q) to the vanishing eigenvalues [9]: in terms of the eigenvectors e Pn i=q+1 eik ein pk = Pn ; k = 1...n . (5) 2 i=q+1 ein The parameter vector m = [p1 . . . pn−1 ]T can be expressed in a vector equation: Pn T i=q+1 ein [ei1 , . . . , ei(n−1) ] Pn m= . (6) 2 i=q+1 ein Thus the solution can be expressed as a linear combination of the ”reduced” ˆ k = √P 1 [ek1 . . . ek(n−1) ]T . The resulting solutions for noreigenvectors: b n−1 2 i=1
eki
mal and full flow on the moving square sequence are given in Fig. 1b and c. However when we are estimating parameters that are of greatly different magnitude some caution is advisable. A scaling procedure will lead to different relative errors of the estimated parameters. Moreover, being a linear combination for q < n−1, the minimum norm solution is dominated by the larger parameters. For pure velocity [6] this does not occur but here it has to be taken into account in further processing steps such as the regularisation procedure described next.
4
Parameter Field Regularisation
We now present a framework that allows to compute dense parameter fields from the previously computed minimum norm solutions. From the structure of ˆk to the the TLS solution given by (6) we can use the reduced eigenvectors b non-vanishing eigenvalues to define a projection matrix which projects onto the subspace determined by the TLS algorithm: P = B q B Tq
where
ˆ1 . . . b ˆq ] . B q = [b
(7)
Each estimated parameter vector m restricts the solution within this subspace. We require the regularised solution r to be close in a least squares sense. A zero weight ω reflects locations where no solution could be found. A conceptually similar data term using the orthogonal projection on a subspace of an ordinary least squares solution based on multiple constraints at the same pixel has been given by [11]. Combining such a data constraint with a simple membrane smoothness constraint, as used by [4], in the the considered area A yields the following minimisation problem (we set ω = ωt in the following): ) Z ( m X 2 2 ω (P r − m) + α (∇ri ) dxdy → min . (8) A
i=1
a
b
c
d
Fig. 2: Moving square optical flow after: a 10, b 100, c 200 and d 500 iterations.
Here α is a regularisation parameter that controls the influence of the smoothness term. Evaluating the Euler-Lagrange equations yields: d d 2ωP (P r − m) − 2α (r x ) + (r y ) = 0 . (9) dx dy The Laplacian ∆r = r xx + r yy can be approximated as ∆r = r¯ − r, where r¯ denotes a local average [12]. Thus we arrive at: (ωP + α1l) r = α¯ r + ωP m .
(10)
This now enables an iterative solution to the minimisation problem. Let A = ωP + α1l, then an update r k+1 from the solution at step k is given by: r k+1 = αA−1 r¯ k + ωA−1 P m .
(11)
Initialisation could be done with zero. However, in order to speed up the convergence we use an interpolation of the available full parameter field by normalised averaging. The matrix A−1 can be found to be: ω −1 −1 A =α 1lm − P . (12) α+ω In order to verify that A−1 A = 1l only the idempotence of the projection matrix is needed. Inserting (12) into (11) yields: r k+1 = r¯ k −
ω ω P r¯ k + Pm . α+ω α+ω
(13)
Using 1l = P + P ⊥ we can separate the two subspaces: pk+1 = P ⊥ r¯ k +
1 P (α¯ r k + ωm) . α+ω
(14)
In the orthogonal subspace, where no initial TLS solution is available, the local average is used. In the TLS subspace the iterative update is given by a weighted mean of the local average and the originally available solution. The weights are determined by the regularisation constant α and the confidence measure ω
a
b
c
d
Fig. 3: IR sequence of the water surface: a original image, b map indicating the type of flow (see Fig. 4), c velocity field and d source term (temperature change) [−0.5, 0.5]K/f rame.
respectively. Notice that the choice of α only regulates the amount of smoothing in the TLS subspace and has no effect in the orthogonal subspace. The effect of this regularisation on the example of Fig. 1 is shown in Fig. 2 (α = 10). Initially the flow field is dominated by the normal flows. With increasing iterations the full flow information spreads out until a dense full flow field is obtained. Here the trace is used to mask the area A. As stated above some of the minimum norm parameters can be unreliable when there is a difference in magnitude between the parameters. This can be accounted for by choosing different confidence values ωi for each parameter ri in (13): ( ωi =
ci ω ω
if q < n − 1 . else
(15)
Where ci is a factor reflecting the confidence that the parameter ri is reliably computed in the minimum norm solutions. In the same way we can also choose a different smoothing strength αi for the different parameters.
5
Experiments
We are now ready to demonstrate the proposed algorithm on a few examples. Heat flux. Using IR cameras to view the water surface gives a direct way to study the heat exchange [8]. Here a linear source term is used to capture the rate of change. An example is given in Fig. 3 where we find that the temperature change (Fig. 3d) does not seem to be correlated to either the original temperature distribution (Fig. 3a) nor the velocity field (Fig. 3c). Diverging tree. Frame 20 of the diverging tree sequence is shown in Fig. 4a [3]. Because we are estimating a second order process we have chosen the following weight and smoothness factors in this case: α1,2 = 1, α3 = 10, c1,2 = 1 and c3 = 0.01. When the divergence of the optical flow field is modeled we obtain a type map as given in Fig. 4b. There is a large region where all parameters can be estimated that corresponds to the tree. Around this area some dependencies occur. The velocity part of the estimated parameter field at these locations
a
b
c
d
e
f
g
h
Fig. 4: Diverging tree sequence: a original image, b map indicating the type of flow, black: tr(J ) < τ1 , white: q = n − 3, light grey: q = n − 2, dark grey: q = n − 1. c velocity corresponding to the first normal flow q = n − 2, d velocity corresponding to the second normal flow q = n − 3. e correct optical flow and f estimated dense optical flow, g correct divergence map and h estimated divergence map [0, 0.04]1/frame.
is shown in Fig. 4c and d. The correct and the regularised velocity are shown in Fig. 4e and f, the only noticeable difference occurs at the image boundaries where the averaging is problematic. The computed divergence (Fig. 4h) is much less accurate, compare to Fig. 4g. However the following numbers indicate that the estimated optical flow does indeed improve when the divergence is modeled, here the angular error Ea is reported [3]: model full flow density only motion p = [u v 1]T 50.2% with div p = [u v d 1]T 73.5%
initial Ea final Ea at 100% 2.4 ± 1.6 1.8 ± 1.3 1.0 ± 0.8 1.0 ± 0.8
When the divergence is included the density of full parameter estimates greatly increases as well as the error drops substantially. The regularisation fills the parameter field without decreasing the optical flow accuracy. Growing maize root. The study of growing plants presents another application where a direct estimation of the divergence (= local growth rate) is useful. As we know that negative growth rates do not make sense we can improve the result by removing the few occurring negative values by setting the confidence to zero. The result obtained by applying our technique to a sequence of a growing maize root is shown in Fig. 5. In the growth map (Fig. 5d) we observe two zones of increased growth: directly at the tip and about 3mm behind the tip.
a
b
c
d
Fig. 5: Growing maize root: a original image, b map indicating the type of flow (see Fig. 4), c movement and d divergence (growth rate) [0, 36]%/h.
6
Conclusion
We introduced a framework that allows to compute dense parameter fields using all the available information from a previous TLS estimation. Besides optical flow this can be used for the extraction of the parameters governing dynamic processes that can be expressed in a linear partial differential equation. The proposed regularisation uses a data constraint that restricts the solution within the subspace determined by the TLS algorithm. Currently a simple membrane smoother is used but we intend to extend this to anisotropic smoothers as well.
Bibliography [1] Haußecker, H., Spies, H.: Motion. In: Handbook of Computer Vision and Applications. Academic Press (1999) [2] Haußecker, H., Garbe, C., Spies, H., J¨ahne, B.: A total least squares for lowlevel analysis of dynamic scenes and processes. In: DAGM, Bonn, Germany (1999) 240–249 [3] Barron, J.L., Fleet, D.J., Beauchemin, S.: Performance of optical flow techniques. International Journal of Computer Vision 12 (1994) 43–77 [4] Horn, B.K.P., Schunk, B.: Determining optical flow. Artificial Intelligence 17 (1981) 185–204 [5] Schn¨ orr, C., Weickert, J.: Variational image motion computation: Theoretical framework, problems and perspectives. In: DAGM, Kiel. Germany (2000) 476–487 [6] Spies, H., J¨ ahne, B., Barron, J.L.: Regularised range flow. In: ECCV, Dublin, Ireland (2000) 785–799 [7] Haußecker, H., Fleet, D.J.: Computing optical flow with physical models of brightness variation. PAMI 23 (2001) 661–673 [8] Garbe, C. S., J¨ ahne, B.: Reliable estimates of the sea surface heat flux from image sequences. In: DAGM, Munich, Germany (2001) 194–201 [9] Van Huffel, S., Vandewalle, J.: The Total Least Squares Problem: Computational Aspects and Analysis. SIAM, Philadelphia (1991) [10] M¨ uhlich, M., Mester, R.: The role of total least squares in motion analysis. In: ECCV, Freiburg, Germany (1998) 305–321 [11] Schn¨ orr, C.: On functionals with greyvalue-controlled smoothness terms for determining optical flow. PAMI 15 (1993) 1074–1079 [12] J¨ ahne, B.: Digital Image Processing. 3 edn. Springer, Germany (1995)
ICG-III: Phytosphere, Research Center J¨ ulich, 52425 J¨ ulich, Germany, [email protected] 2 Interdisciplinary Center for Scientific Computing, University of Heidelberg, INF 368, 69120 Heidelberg, Germany, [email protected]
Abstract. A method for the interpolation of parameter fields estimated by total least squares is presented. This is applied to the study of dynamic processes where the motion and further values such as divergence or brightness changes are parameterised in a partial differential equation. For the regularisation we introduce a constraint that restricts the solution only in the subspace determined by the total least squares procedure. The performance is illustrated on both synthetic and real test data.
1
Introduction
The parameters governing the changes in an image sequence can be estimated locally by means of a total least squares technique [1; 2]. For the simple case of movement alone, i.e. optical flow, it is known that in areas with one dominant gradient direction only the component of the velocity in that direction can be computed. This is the well known aperture problem [3]. When there are more parameters to be estimated there might be even more such linear dependencies between them. We show how such cases can be detected and how an appropriate minimum norm solution can be obtained. In many applications we are nevertheless interested in dense parameter fields and thus there is the need for an interpolation procedure. Clearly it is desirable to utilise the available information from the minimum norm solutions. Towards this end we present a regularisation framework that restricts the interpolated solution only in the subspace in which a solution is available. This is in contrast to variational methods that directly seek to minimise the constraint over the entire data set [4; 5]. The main idea of the presented regularisation has been introduced for the special case of 3D velocity (range flow) estimation in [6]. Here this is extended and applied to more general dynamic processes.
2
Dynamic Processes in Image Sequences
In scientific image processing one often encounters a linear equation for a set of model parameters m at each data point (pixel) of type: aT m = b; a, m ∈ IRm .
An ordinary least squares estimate implies that only b is erroneous while a is exact. This assumption certainly does not hold for the cases discussed below where a contains derivatives. Therefore we introduce a new parameter vector p of dimension n = m + 1 and rewrite the constraint as: dT p = 0 with d = (aT , −b)T
and p = (mT , 1)T .
(1)
This concept is very general in the sense that the parameters of any dynamic process that can be modeled by a linear partial differential equation falls into this scheme [2; 7]. Prominent physical examples are source terms, relaxation and diffusion processes. If there is no such model available a Taylor series expansion up to the desired order may be used [7]. A few examples that are used in the later experiments (Sect. 5) are: Optical flow. The assumption that all temporal variation in the image brightdy dx ness g(x, y, t) is caused by movements implies [4]: dg dt = gx dt + gy dt + gt = 0. With indices denoting partial derivatives. This brightness change constraint equation forms the basis of many algorithms to compute the optical flow dy T [u, v]T = [ dx dt , dt ] [3; 1]. This is readily expressed in the form of (1) using T d = [gx gy gt ] and p = [u v 1]T . Source terms. Source terms cause first-order temporal changes in the image seT quence proportional to a source strength q: ( dg dt = q) yielding d = [gx gy 1 gt ] T and p = [u v q 1] . Source terms occur for example in infrared images if an object is heated or cooled or if the global illumination of a scene changes [8]. The same equation is obtained in the study of the three-dimensional movement of surfaces from sequences of depth maps: zx u + zy v + zt = w. Here [u v w]T denotes the 3D-movement also called range flow [6]. Divergent velocity. The proposed concept can also model spatial variations in the velocity field, e.g. an affine velocity field [1]. If we are interested in a particular aspect of the flow field, such as its divergence or rotation, we can only include this in our estimation. For divergence d this reads as: d = [gx gy (xgx +ygy ) gt ]T ; p = [u v d 1]T . Here x, y denote local coordinates with respect to the computed location. Examples include camera movement along the optical axis or local growth rates in biological experiments.
3
TLS Estimation
The presented estimation framework extends previous work on total least squares optical flow estimation [2]. For an estimation we pool the constraints (1) over a local neighbourhood and assume the parameters to be constant within: Z ∞ ˆ = arg min ˆ )2 dx0 dt0 subject to p ˆT p ˆ = 1 . (2) p w(x − x0 , t − t0 )(dT p −∞
Here w(x − x0 , t − t0 ) is a weighting function that defines the spatio-temporal neighbourhood where the parameters are estimated. In order to avoid the trivial
a
b
c
d
Fig. 1: Moving square example: a image data, b type of flow (black: trace to small, light grey: linear dependency and dark grey: full estimate possible), c normal flow and d full flow.
ˆ to be normalised. Minimisation leads solution we require the solution vector p to an eigenvalue equation of the extended structure tensor J : Z ∞ ˆ = λˆ Jp p with J (x) = w (x − x0 , t − t0 ) (d dT ) dx0 dt0 . (3) −∞
It follows that the solution to (2) is given by the nullspace of J [9]. Section 3.2 describes how to recover the parameters. For an unbiased and optimal TLS estimate the errors in all data terms need to be independent random variables with zero mean and equal variance [9; 10]. This implies a scaling with respect to the absolute error values. 3.1
Type and Confidence Values
The spectrum of J can be used to detect linear dependencies in the data. In the following we assume the eigenvalues of J to be sorted in descending order: λ1 ≥ λ2 ≥ . . . ≥ λn−1 ≥ λn . The trace of J is used to determine if there is enough variation (texture) in the data for the computation to make sense. Thus we only proceed where tr(J ) > τ1 . The smallest eigenvalue directly yields the residual of the parameter fit. Thus we reject unreliable estimates if λn > τ2 , with a threshold τ2 corresponding to the noise level. This occurs at discontinuities where a single set of parameters is not sufficient or when no coherent motion is present. A small λn does not ensure that all the parameters can be estimated independently from each other. More than one eigenvalue may be close to zero (< τ2 ) and we can no longer uniquely pick a solution. Any vector in the nullspace of J is a possible solution, see Sect. 3.2. A type measure that indicates how well the dimensionality of the nullspace of J has been determined can be obtained by examining how much the first relevant eigenvalue λq is above the threshold: 2 λ q − τ2 ωt = , λ1 ≥ . . . ≥ λq > λq+1 ≈ . . . ≈ λn < τ2 . (4) λq For a simple optical flow example of a moving square (Fig. 1) a map indicating which type of flow can be computed is given in Fig. 1b.
3.2
Minimum Norm Solutions
As stated above any vector in the nullspace of J is a possible solution. A sensible choice is that resulting in an estimate with minimal norm. This solution is given ˆi (i > q) to the vanishing eigenvalues [9]: in terms of the eigenvectors e Pn i=q+1 eik ein pk = Pn ; k = 1...n . (5) 2 i=q+1 ein The parameter vector m = [p1 . . . pn−1 ]T can be expressed in a vector equation: Pn T i=q+1 ein [ei1 , . . . , ei(n−1) ] Pn m= . (6) 2 i=q+1 ein Thus the solution can be expressed as a linear combination of the ”reduced” ˆ k = √P 1 [ek1 . . . ek(n−1) ]T . The resulting solutions for noreigenvectors: b n−1 2 i=1
eki
mal and full flow on the moving square sequence are given in Fig. 1b and c. However when we are estimating parameters that are of greatly different magnitude some caution is advisable. A scaling procedure will lead to different relative errors of the estimated parameters. Moreover, being a linear combination for q < n−1, the minimum norm solution is dominated by the larger parameters. For pure velocity [6] this does not occur but here it has to be taken into account in further processing steps such as the regularisation procedure described next.
4
Parameter Field Regularisation
We now present a framework that allows to compute dense parameter fields from the previously computed minimum norm solutions. From the structure of ˆk to the the TLS solution given by (6) we can use the reduced eigenvectors b non-vanishing eigenvalues to define a projection matrix which projects onto the subspace determined by the TLS algorithm: P = B q B Tq
where
ˆ1 . . . b ˆq ] . B q = [b
(7)
Each estimated parameter vector m restricts the solution within this subspace. We require the regularised solution r to be close in a least squares sense. A zero weight ω reflects locations where no solution could be found. A conceptually similar data term using the orthogonal projection on a subspace of an ordinary least squares solution based on multiple constraints at the same pixel has been given by [11]. Combining such a data constraint with a simple membrane smoothness constraint, as used by [4], in the the considered area A yields the following minimisation problem (we set ω = ωt in the following): ) Z ( m X 2 2 ω (P r − m) + α (∇ri ) dxdy → min . (8) A
i=1
a
b
c
d
Fig. 2: Moving square optical flow after: a 10, b 100, c 200 and d 500 iterations.
Here α is a regularisation parameter that controls the influence of the smoothness term. Evaluating the Euler-Lagrange equations yields: d d 2ωP (P r − m) − 2α (r x ) + (r y ) = 0 . (9) dx dy The Laplacian ∆r = r xx + r yy can be approximated as ∆r = r¯ − r, where r¯ denotes a local average [12]. Thus we arrive at: (ωP + α1l) r = α¯ r + ωP m .
(10)
This now enables an iterative solution to the minimisation problem. Let A = ωP + α1l, then an update r k+1 from the solution at step k is given by: r k+1 = αA−1 r¯ k + ωA−1 P m .
(11)
Initialisation could be done with zero. However, in order to speed up the convergence we use an interpolation of the available full parameter field by normalised averaging. The matrix A−1 can be found to be: ω −1 −1 A =α 1lm − P . (12) α+ω In order to verify that A−1 A = 1l only the idempotence of the projection matrix is needed. Inserting (12) into (11) yields: r k+1 = r¯ k −
ω ω P r¯ k + Pm . α+ω α+ω
(13)
Using 1l = P + P ⊥ we can separate the two subspaces: pk+1 = P ⊥ r¯ k +
1 P (α¯ r k + ωm) . α+ω
(14)
In the orthogonal subspace, where no initial TLS solution is available, the local average is used. In the TLS subspace the iterative update is given by a weighted mean of the local average and the originally available solution. The weights are determined by the regularisation constant α and the confidence measure ω
a
b
c
d
Fig. 3: IR sequence of the water surface: a original image, b map indicating the type of flow (see Fig. 4), c velocity field and d source term (temperature change) [−0.5, 0.5]K/f rame.
respectively. Notice that the choice of α only regulates the amount of smoothing in the TLS subspace and has no effect in the orthogonal subspace. The effect of this regularisation on the example of Fig. 1 is shown in Fig. 2 (α = 10). Initially the flow field is dominated by the normal flows. With increasing iterations the full flow information spreads out until a dense full flow field is obtained. Here the trace is used to mask the area A. As stated above some of the minimum norm parameters can be unreliable when there is a difference in magnitude between the parameters. This can be accounted for by choosing different confidence values ωi for each parameter ri in (13): ( ωi =
ci ω ω
if q < n − 1 . else
(15)
Where ci is a factor reflecting the confidence that the parameter ri is reliably computed in the minimum norm solutions. In the same way we can also choose a different smoothing strength αi for the different parameters.
5
Experiments
We are now ready to demonstrate the proposed algorithm on a few examples. Heat flux. Using IR cameras to view the water surface gives a direct way to study the heat exchange [8]. Here a linear source term is used to capture the rate of change. An example is given in Fig. 3 where we find that the temperature change (Fig. 3d) does not seem to be correlated to either the original temperature distribution (Fig. 3a) nor the velocity field (Fig. 3c). Diverging tree. Frame 20 of the diverging tree sequence is shown in Fig. 4a [3]. Because we are estimating a second order process we have chosen the following weight and smoothness factors in this case: α1,2 = 1, α3 = 10, c1,2 = 1 and c3 = 0.01. When the divergence of the optical flow field is modeled we obtain a type map as given in Fig. 4b. There is a large region where all parameters can be estimated that corresponds to the tree. Around this area some dependencies occur. The velocity part of the estimated parameter field at these locations
a
b
c
d
e
f
g
h
Fig. 4: Diverging tree sequence: a original image, b map indicating the type of flow, black: tr(J ) < τ1 , white: q = n − 3, light grey: q = n − 2, dark grey: q = n − 1. c velocity corresponding to the first normal flow q = n − 2, d velocity corresponding to the second normal flow q = n − 3. e correct optical flow and f estimated dense optical flow, g correct divergence map and h estimated divergence map [0, 0.04]1/frame.
is shown in Fig. 4c and d. The correct and the regularised velocity are shown in Fig. 4e and f, the only noticeable difference occurs at the image boundaries where the averaging is problematic. The computed divergence (Fig. 4h) is much less accurate, compare to Fig. 4g. However the following numbers indicate that the estimated optical flow does indeed improve when the divergence is modeled, here the angular error Ea is reported [3]: model full flow density only motion p = [u v 1]T 50.2% with div p = [u v d 1]T 73.5%
initial Ea final Ea at 100% 2.4 ± 1.6 1.8 ± 1.3 1.0 ± 0.8 1.0 ± 0.8
When the divergence is included the density of full parameter estimates greatly increases as well as the error drops substantially. The regularisation fills the parameter field without decreasing the optical flow accuracy. Growing maize root. The study of growing plants presents another application where a direct estimation of the divergence (= local growth rate) is useful. As we know that negative growth rates do not make sense we can improve the result by removing the few occurring negative values by setting the confidence to zero. The result obtained by applying our technique to a sequence of a growing maize root is shown in Fig. 5. In the growth map (Fig. 5d) we observe two zones of increased growth: directly at the tip and about 3mm behind the tip.
a
b
c
d
Fig. 5: Growing maize root: a original image, b map indicating the type of flow (see Fig. 4), c movement and d divergence (growth rate) [0, 36]%/h.
6
Conclusion
We introduced a framework that allows to compute dense parameter fields using all the available information from a previous TLS estimation. Besides optical flow this can be used for the extraction of the parameters governing dynamic processes that can be expressed in a linear partial differential equation. The proposed regularisation uses a data constraint that restricts the solution within the subspace determined by the TLS algorithm. Currently a simple membrane smoother is used but we intend to extend this to anisotropic smoothers as well.
Bibliography [1] Haußecker, H., Spies, H.: Motion. In: Handbook of Computer Vision and Applications. Academic Press (1999) [2] Haußecker, H., Garbe, C., Spies, H., J¨ahne, B.: A total least squares for lowlevel analysis of dynamic scenes and processes. In: DAGM, Bonn, Germany (1999) 240–249 [3] Barron, J.L., Fleet, D.J., Beauchemin, S.: Performance of optical flow techniques. International Journal of Computer Vision 12 (1994) 43–77 [4] Horn, B.K.P., Schunk, B.: Determining optical flow. Artificial Intelligence 17 (1981) 185–204 [5] Schn¨ orr, C., Weickert, J.: Variational image motion computation: Theoretical framework, problems and perspectives. In: DAGM, Kiel. Germany (2000) 476–487 [6] Spies, H., J¨ ahne, B., Barron, J.L.: Regularised range flow. In: ECCV, Dublin, Ireland (2000) 785–799 [7] Haußecker, H., Fleet, D.J.: Computing optical flow with physical models of brightness variation. PAMI 23 (2001) 661–673 [8] Garbe, C. S., J¨ ahne, B.: Reliable estimates of the sea surface heat flux from image sequences. In: DAGM, Munich, Germany (2001) 194–201 [9] Van Huffel, S., Vandewalle, J.: The Total Least Squares Problem: Computational Aspects and Analysis. SIAM, Philadelphia (1991) [10] M¨ uhlich, M., Mester, R.: The role of total least squares in motion analysis. In: ECCV, Freiburg, Germany (1998) 305–321 [11] Schn¨ orr, C.: On functionals with greyvalue-controlled smoothness terms for determining optical flow. PAMI 15 (1993) 1074–1079 [12] J¨ ahne, B.: Digital Image Processing. 3 edn. Springer, Germany (1995)
