Motion analysis from first-order properties of optical flow

CVGIP:

IMAGE

UNDERSTANDING

Vol. 56, No. 1, July, pp. 90-107, 1992

Motion Analysis from First-Order Properties of Optical Flow MARCO CAMPANI

AND ALESSANDRO

VERRI

Dipartimento di Fisica dell’Universitd di Genova, Via Dodecaneso 33, 16146 Genoa, Italy Received April 10, 1991; accepted March 18, 1992

In this paper the first-order spatial propertiesof optical flow, such as singular points and elementarymotions, are usedto describe and analyze moving images. First, it is shown that the spatial structure of the viewed motion of rigid objects,with the exceptionof surfaceand motion boundaries,can usually be approximated over rather large regions of the image plane by a linear vector field. Second,a methodis proposedin which optical flow is computedas a piecewiselinear vector field. The obtained first-order propertiesof optical flow are then usedto devisemethodsfor distinguishing betweendifferent kinds of motion, extracting qualitative information about shape, and segmentingthe viewed image into the different moving objects.Experimental resultson severalsequencesof real imagesare presentedand discussed.Sincethe methodswhich are proposedrely upon patchwise and not pointwise flow estimates,the obtained results are usually very goodand insensitiveto noise.It is concludedthat the first-orderpropertiesof optical flow are very helpful for the understanding of ViSud motion. 0 1992 Academic PWSS, h.

1. INTRODUCTION In recent years the analysis of moving images has received increasing attention from the computer vision community. A number of feature-based methods have been proposed for the reconstruction of the three-dimensional (3D) motion and structure of the viewed scene from time-varying images [l-5]. In essence, these methods view the apparent motion of each feature as the perspective projection of the corresponding motion in space and then solve the equations for the 3D motion and structure by means of additional assumptions on the rigidity of the observed motion and the geometry of the viewed structure. A common property of these methods is that the estimates of image motion of sparse and well defined features can be quite accurate. However, feature-based analysis does not appear to capture the underlying coherence of the observed motion, nor to cope efficiently with scenarios in which different objects may be moving independently. Instead, a somewhat alternative point of view consists of finding strategies that look at the apparent motion of every point of the image plane [6-131. The visual infor-

mation can then be thought of as encoded in a dense twodimensional (2D) vector field, called optical flow [14], which describes the motion of the time-varying image brightness pattern. The main advantages of this second type of strategy are that optical flow enhances the spatial coherence underlying the motion of nearby image points and can be computed independently of the complexity of the viewed scene. Ambiguity in the definition of optical flow and high sensitivity to noise of the estimates of image brightness derivatives are usually reported as the major shortcomings. In this paper it is shown that first-order spatial properties of optical flow can be effectively employed for the solution of problems such as distinguishing between different kind of motions, extracting qualitative information about shape, and segmenting the viewed image in the different moving objects. The key point is that, over rather large patches of the image plane, the optical flow of rigid objects can usually be approximated by a linear vector field. Therefore, first-order spatial differential properties of optical flow, such as singular points (or the points where optical flow vanishes) and elementary motions (such as translation, expansion, and rotation), are essentially sufficient to describe image motion. Moreover, optical-flow-based methods can rely upon patchwise, instead of pointwise, flow estimates. This overcomes the problems both of ambiguity in the definition of optical flow and of sensitivity of image brightness derivatives to noise. The paper is organized as follows. Section 2 analyzes the approximation of optical flow in terms of its firstorder spatial properties. These properties, which are the singular points and elementary motions, are then briefly discussed in Section 3. Section 4 describes an algorithm for the computation of optical flow which takes advantage of the simple piecewise linear structure of optical flow. The extraction of qualitative information from firstorder spatial properties of optical flow of real images is described in Sections 5, 6, and 7. Section 5 is devoted to 3D motion estimation, Section 6 to surface reconstruction from translation, and Section 7 to image segmentation. Finally, Section 8 summarizes the obtained results.

90 1049-9660192$5.00 Copyright 0 1992 by Academic Press, Inc. A11rights of reproduction in any form reserved

MOTION ANALYSIS

91

FROM OPTICAL FLOW

2. SPATIAL STRUCTURE OF OPTICAL FLOW -

The aim of this section is to substantiate the empirical observation that the spatial structure of a typical optical flow, apart from the discontinuities which generate from surface and motion boundaries, appears to be almost always “fairly simple.” The optical flow of moving planes is first studied. The analysis is then extended to nonplanar surfaces. 2.1.

pfX x3

’

where Sis the focal length. Note that x3 = f, for every x. If the point X is moving with velocity V, then the velocity v of the corresponding point x can be obtained by differentiating Eq. (1) with respect to time; that is, e3 X (V x X) x: .

The vector field v in Eq. (2) is called a motionfield. Since u3 is always zero, we write v = (ui , uz). The term optical flow usually refers to an estimate of the motion field which can be obtained from some suitable combination of the spatial and temporal changes of the image brightness pattern E = E(x, t), where t is the time. Thus, whereas the motion field is a mathematical concept which is unambiguously defined, the notion of optical flow may vary depending on the assumptions which are used for the computation of the optical flow itself. In the following, however, for the sake of simplicity, the term optical flow will also be used to denote the motion field. The context will always make it easy to distinguish between these two somewhat different concepts. 2.2.

CZl3XlX2 -

a12fx1

+

a3d2 U23.X:

-

(3) +

(U33

-

U22)fX2

a32f2,

where the aij, i, j = I, 2, 3, are given by (I) . ej X ei ai.i=fs (g(t)

UiTj),

(4)

Preliminaries

Let us first define the main concepts and notations which will be used throughout the paper. In what follows if w is an arbitrary vector, we let w;, i = 1, 2, 3, be the components of w in S, a 3D system of orthogonal coordinates, and e;, i = 1, 2, 3, the three unit vectors parallel to the orthogonal axes. The image plane is orthogonal to e3 (that is, the optical axis is parallel to e3), while the focus of projection lies at the origin 0 of S. Now let X be a point in 3D space and x its perspective projection onto the image plane. We then have

v=f

U2 =

a2dx2

The Optical Flow of Planar Surfaces Revisited

Let us now consider the optical flow of a planar surface. If the 3D velocity of a point X is written as V = w x X + T, where 6~and T are the angular and translational velocity respectively, the optical flow components u, and v2 can be written as

where u is the unit normal vector of the planar moving surface which, at each time t, satisfies the equation g(t) = u . x.1

The optical flow components of Eqs. (3) are usually regarded as quadratic polynomials in the x1 and x2 coordinates. However, the same components can also be regarded as homogeneous quadratic polynomials in the xi, x2, and x3 coordinates (as x3 = f). Since f is usually larger than L, where L is the linear size of the viewing sensor, f is also larger than both xi and x2. This last observation is now used to show that the optical flow of a planar surface is usually well approximated, over rather large areas of the image plane, by a linear vector field. Let us divide the optical flow of a moving planar surface into a number of squared patches of linear size equal to s and compute the best linear approximation of the optical flow in each patch. Apart from special values of the motion and structure parameters for which some of the coefficients av may vanish, it can be expected that, for sufficiently small values of s, the obtained piecewise linear approximation will be accurate. Since Eq. (3) is homogeneous in the image coordinates xl , x2, and f, it is easy to show that the accuracy of the approximation actually depends on the ratio r = s/f. Since fis larger than both xi and x2 for each point x of the image plane, this ratio is small even for relatively large values of s. Let us further discuss this point by using an interesting example. Figure 1A shows the optical fow of a planar surface which rotates around a fixed vertical axis. The focal length f is twice as long as the linear size L over which optical flow has been computed. The optical flow of Fig. IA is highly nonlinear due to the presence of three singular points within the field of view. Since the optical flow of a planar surface can have at most three singular points Ill], the complexity of the spatial structure of the flow of Fig. 1A is somehow maximal. Figure IB shows the piecewise linear optical flow, which can be obtained by dividing the optical flow of Fig. IA into 16 squares of linear size s = L/4 and computing the linear vector field ’ Note that from the definition of the a’J in Eq. (4) the classical twofold ambiguity of the optical flow of planar surfaces [15-1X] easily follows for the a,, and thus Eqs. (3) are left unchanged if u + TIT, T -+ Tu, and

w+ w + u x T/g(r).

92

CAMPANI AND VERRI

60

0 20

30

40

SO 60

70

80

90108 0

FIG. 1. (A) Subsampled optical flow of a planar surface rotating around a fixed axis. The focal length f equals 1000, while the motion and structure parameters are w = (0, 1, 0), T = (- lOOf, 0, 0), II = (.07, -. 11, .99), and g = 98f. (B) The optical flow which is obtained by dividing the field of view (512 x 512 pixels) into 16 squared patches and computing the linear field which best approximates the optical flow of (A) in each patch. (C) and (D) Histograms of the directions and speeds, respectively, of the optical flows of (A) (white) and (B) (textured columns). The histograms are normalized to the maximum value.

which best approximates the original flow in each square through standard least-mean-square techniques. The similarity between the two flows is remarkable (the ratio Y = s/f equals l/8). The two histograms of Figs. 1C and D, respectively, compare the direction and speed distributions of the optical flows of Figs. 1A and B. It is evident that, even at this finer detail, the difference between the quadratic and the piecewise linear vector field is very small. Extensive experimentation on examples of comparable complexity gave similar results. It can thus be concluded that the optical flow of a planar surface can usually be well approximated, over rather large areas of the image plane, by a linear vector field. 2.3.

General Case

Let us now extend the previous analysis to the case of nonplanar surfaces by discussing a simple example. Fig-

ure 2A shows the optical flow of a sphere which translates toward the image plane. The focal length f is again twice as long as the linear size L of the optical flow. The optical flow of Fig. 2A is highly nonlinear and cannot be written as a polynomial of any finite order in the image coordinates. As in the previous example, Fig. 2B shows the piecewise linear optical flow which can be obtained by dividing the optical flow of Fig. 2A into 16 squares of linear size s = L/4 and computing the linear vector field which best approximates the original flow in each square through standard least-mean-square techniques. The two flows are very similar, as can be inferred from the histograms of Figs. 2C and D in which the direction and speed distributions of the optical flows of Figs. 2A and B are respectively compared. Two reasons account for this fact. First, as in the previous example, the accuracy of the linear approximation increases for small values of the ratio r (although in this case the optical flow components

MOTION ANALYSIS

100

100

80

80

60

60

40

40

20

20

0

0

1

1 90

1

1 180

1

1 270

1

93

FROM OPTICAL FLOW

1

0

0

10

20

30

40

50

60

70

80

90100

FIG. 2. (A) Subsampled optical flow of a sphere translating along the optical axis. The focal length f equals 1000, the motion parameters are w = (0, 0,O) and T = (-lOof, 0, 0), while the sphere radius is lflfand the distance between the sphere center and the image plane is IOllf. (B) The optical flow which is obtained by dividing the field of view (512 X 512 pixels) into 16 squared patches and computing the linear field which best approximates the optical flow of (A) in each patch. (C) and (D) Histograms of the directions and speeds, respectively. of the optical flows of (A) (white) and (B) (textured columns). The histograms are normalized to the maximum value.

are not homogeneous polynomials in the image coordinates). Second, and most importantly, large deviations from linearity can be expected where large changes in the spatial properties of the viewed surface occur over relatively small areas of the image plane. These changes, however, usually correspond to image boundaries where optical flow is more likely to be discontinuous. To summarize, it appears that, due to the relatively small aperture of a typical imaging device with respect to the focal length, the optical flow of rigid moving objects, with the exception of image boundaries, can usually be well approximated, over rather large areas of the image plane, by a linear vector field. In Section 4 this observation will be used as a constraint for the computation of optical flow. 3.

FIRST-ORDER SPATIAL PROPERTIES OPTICAL FLOW

OF

Since the spatial structure of optical flow is usually well described by a piecewise linear vector field, the first-

order spatial properties of optical flow, such as singular points and elementary motions, are likely to be relevant features for both the computation and the analysis of optical flow. This section briefly reviews the main definitions and properties of the singular points and elementary motions of optical flow. 3.1.

Singular Points

The linear part of optical flow in the neighborhood of a point x can be described by the 2 x 2 matrix M whose entries are

where the partial derivatives are computed at x. If x is a singular point, that is, v(x) = 0, the matrix M of Eq. (5) can be useful to distinguish between different kind of motions and to estimate motion parameters. If x is not a singular point, the matrix M can be used to understand

94

CAMPANI

qualitative properties of the viewed motion and scene, and to segment the image into the different moving objects. Let us first consider the case in which x is a singular point. The qualitative nature of 3D motion (like translation or rotation) can be described by looking at the temporal evolution of the spatial structure of optical flow in the neighborhood of the singular points [ 1I]. We first look at the simple but interesting case of translational motion. 3.1.1.

Translation

Let S be a surface which is translating with velocity V = T. The resulting optical flow v has at most one singular point p such that [1 11 p=f.$

(6)

From Eq. (6), it can easily be seen that 1. The spatial location of p on the image plane does not change over time. 2. The heading direction is related to the location of p over the image plane. 3. The matrix M at p is a multiple of the identity matrix; that is, M = Al. In the language of the theory of dynamical systems, p is always a focus [ 191with the only eigenvalue X of M given by

AND VERRI

angular velocity o can be written in terms of the complex conjugate eigenvalues h and h of the matrix M computed at p as 62 = Ah. In all the other cases, the spatial structure of the flow in the neighborhood of an immobile point p changes over time depending on the angle 8 between the rotation axis and the unit normal to the surface at p and on the relative position of the rotating surface with respect to the viewing point (see [ 1I]). 3.1.3.

General Motion

The case of general motion is more complicated. Let us restrict the present analysis of general motion to a simple but interesting case in which the viewed camera is moving on a flat surface S. In this case, which can be termed passive navigation, the rotation axis is orthogonal to both the surface S and the translational component of motion, and useful information can again be obtained from the singular points of optical flow. First, the 3D motion is instantaneously indistinguishable from a rotation around a fixed axis orthogonal to S. Therefore, the singular points must lie on the straight line of the image plane which is parallel to S and contains the optical axis [ 111. At a singular point p we can write

(7)

A=-$$,

7

where r is the time-to-collison. From these properties it is evident that useful information on translational motion can easily be obtained from the singular point of optical flow. 3.1.2. Rotation In the case of a rotating surface it is useful to distinguish between two different kinds of singular points. Let p be the singular point perspective projection of a point P in the 3D space and V(P) the velocity of P. A point p can be a singular point either because P lies on the rotation axis (that is, v(p) = 0 and V(P) = 0), or because V(P) lies on the straight line which goes through P and the center of projection (that is, v(p) = 0 but V(P) # 0). To avoid ambiguities a singular point of the first kind is named an immobile point. Let use briefly discuss a particular case in which p is an immobile point (a more general analysis of the singular and immobile points of rotation can be found in [l 11). If the rotation axis is orthogonal to the plane tangential to the surface at P, p is always a center [19] and the flow in the neighborhood of p is tangential to closed orbits. The

where Det M and Tr M are the determinant and trace of the matrix M, respectively, II is the unit vector normal to the moving surface at P, o is the angular velocity, and V = V(P). From Eq. (7) it follows that the eigenvalues of the matrix M are always real and thus that the singular points cannot be spirals [19]. If some a priori information on the orientation of the viewed surface is available, it can be shown that the motion trajectory and parameters can be accurately estimated from the spatial structure and time evolution of the singular points [20]. Let us now discuss the case in which the matrix M of Eq. (5) is computed at x, where x is not a singular point. 3.2.

Elementary Motions

According to a classical theorem of the theory of deformable bodies [21], the spatial structure of a vector field over a sufficiently small patch can be described as the sum of a rigid translation, a uniform expansion, a pure rotation, and two components of shear. The relevance of this theorem to the analysis of visual motion depends upon the simple piecewise linear spatial structure of optical flow. In Section 2 it has been shown that the optical flow v, at each point x of a neighborhood (not necessarily

MOTION

ANALYSIS

FROM OPTICAL

VE * v + E, = 0,

very small) of a point x0 of the image plane, can usually be approximated as a linear vector field v = v. + M(x - x0),

(8)

where v. is the optical flow at x0 and the matrix M is computed at x0. It is clear that in Eq. (8) the rigid translation is given by the term vo. The matrix M can be written as

95

FLOW

(9)

where V is the two-dimensional spatial gradient and E, the partial temporal derivative of the image brightness E = E(x, t) at the location x on the image plane and at the time t. It has been proved [25] that, with the exception of poorly textured regions and regions with strong highlights, Eq. (9) can be regarded as approximately correct. Section 2 has shown that over rather large areas of the image plane optical flow can be well approximated by the linear vector field of Eq. (8). Thus, by combining Eqs. (8) and (9) it is easy to obtain

where VE . v. + VE * M(x - x0) = -Et.

and CY= (MI1 + M&/2, p = (MI2 - Md2, y = (MI1 M&/2, and 6 = (MI1 + M2J/2. The matrices Ii, i = 1, . . . ) 4, correspond to the elementary motions of uniform expansion, pure rotation, and two components of shear, respectively. Since I, and 12are invariant for arbitrary orthogonal transformation, the quantities cx and p are independent of the system of coordinates. Consequently, the amounts of uniform expansion and pure rotation are intrinsic properties of optical flow [9]. In contrast, since I3 and I4 are not invariant, y and 6 do not describe intrinsic properties of optical flow. An invariant measure of the amount of shear can be given by 7 = V$V3 (since ff2 + p2 + y2 + a2 is invariant). In the next sections both the singular points and the elementary motions are used to identify different kind of motions, estimate motion parameters, extract qualitative information on shape from translational motion, and segment the viewed image in the different moving objects from real images.

(10)

Equation (10) can be thought of as a linear equation for v. and M (that is, for six unknowns). Note that while VE and Et, the equation coefficients, are computed at x, the unknowns, v. and M, refer to the point x0. An algorithm for the computation of optical flow can thus be implemented as follows [22]: Cover the image plane with n x II equally spaced and possibly overlapping squared patches of N x N pixels. Let x6 be the central point of the ith patch, and vi and Mi the optical flow and the matrix of the flow first-order spatial derivatives at x6, i = 1, . . . , II x n. Compute the first-order spatial and temporal derivatives of the image brightness pattern at every point x. Solve, by means of standard least-mean-square techniques, n X n overconstrained systems of 6 x N equations such as Eq. (10) for vb and Mi at the ith patch, i = I, . . . ) n x n; l

l

l

It is apparent that, in order to enforce spatial coherence in the estimation of optical flow, N, the patch linear size, should be as large as possible. On the other hand, n, the patch linear number, controls the redundancy in the estimation of optical flow, redundancy that may be useful to improve accuracy and achieve robustness. Experiments on sequences of synthetic and real images have shown that on 256 x 256 images accurate and robust 4. A METHOD FOR THE COMPUTATION OF estimation can be obtained when N varies between 10 OPTICAL FLOW and 80 pixels, and n between 10 and 70. Note that the In this section a method which has been recently pro- range over which N varies is in very good agreement with posed for the computation of optical flow [22] is de- the results of the analysis of Section 2 on the linear size of scribed. The method, which is demonstrated on se- the patch within which optical flow can be approximated quences of noise corrupted synthetic images, exploits the by a linear vector field. piecewise linear spatial structure of the optical flow discussed in Section 2. 4.2. Experimental Tests 4.1.

Description of the Algorithm

Let us assume that at each location of the image plane the image brightness E changes according to the equation

W, 241

Let us now present some results which have been obtained by applying the described algorithm to sequences of synthetic images. The major aim of the experiments was to test the accuracy and sensitivity to noise of the proposed method.

96

CAMPANI AND VERRI

c

D

8000

,,

,.,,,,,

,,,,,

,

,

,,

,

,,,

,

,

,,

,

r

8000

,(

1000

7000

6000

6000

5000

5000

4000

4000

3000

3000

2000

2000

1000

1000 0

0 0

32

64

96

!8

160

192

224

; 15c

0

A 32

64

96

1

192

224

256

FIG. 3. (A) The “plaid” pattern used in the experimental tests. The image brightness E at the point (xi, x2) is given by the superimposition of two sinusoidal waves; that is, E(x, ,x2) = 128 + 40 sin(w, x1) + 40 cos(wlxJ with WI = w2 = .4. (B) The same plaid pattern as (A) corrupted by white Gaussian noise with standard deviation on = 5. (C) and (D) Gray-value distributions of the patterns of (A) and (B), respectively.

Figures 3A and B show the “plaid” patterns (two superimposed horizontal and vertical sinusoidal waves of variable frequency and amplitude) which have been used to evaluate the performances of the algorithm. The plaid pattern of Fig. 3A is noise-free while the pattern of Fig. 3B is corrupted by a large amount of noise (the distributions of the grey values of the two images of Figs. 3A and B are shown in Figs. 3C and D, respectively). In the first set of experiments we have considered pure translation parallel to the image plane. Table 1 reports the results of experiments in which the plaid pattern was translating in different directions. From Table 1 it can easily be seen that the quantitative agreement between the true and the estimated displacement is always very good and almost insensitive to the added noise. In the second set of experiments we have considered translations of the plaid pattern in the direction orthogonal to the image plane. Figure 4A shows the optical flow obtained while the noise-corrupted pattern was translating toward the image plane. The histogram in Fig. 4B reproduces the average amount of the estimated expansion, rotation, and two components of shear. It is evident

that both the optical flow and the elementary motions have been recovered correctly despite the rather large amount of added noise. This is also confirmed by the histograms of Figs. 4C and D, where the direction and speed distributions of the optical flow of Fig. 4A and of the noise-free corresponding flow are shown. In the third set of experiments we have considered rotation on the image plane. Figure .5A reproduces the

TABLE

1

Estimatesof the AverageOptical Flow Components~1and v2 for Three Cases of Pure Translation

“, = 10.0 cl>= 0.0 u, = 0.0 v2 = 10.0 u, = 10.0 u* = 10.0

01

UUI

u2

VU!

9.98 0.02 10.01

0.02 0.06 0.11

0.01 9.99 10.00

0.07 0.02 0.11

Note. The true values are shown in the leftmost column. vu, and oUz are the mean square errors of u, and u2respectively; (T, = 5 is the noise standard deviation; and cr, = 1 is the standard deviation of the Gaussian filter which has been used to compute image derivatives.

MOTION ANALYSIS

97

FROM OPTICAL FLOW

80 60

0

a

P

Y

6

80 60

0

10

30 40 50

20

60

70

80

90 10

FIG. 4. (A) The optical flow obtained through the proposed method while the noisy pattern of Fig. 3 was translating toward the image plane (N = 41 and n = 42). (B) Average estimates of the amount of expansion, rotation, and two components of shear. (C) Direction distribution and (D) speed distribution for the optical flow obtained from the noise-free sequence (white columns) and the optical flow of (A) (textured columns). The histograms are normalized to the maximum value.

optical flow obtained by using the proposed method while the noise-corrupted pattern was rotating around a fixed point at the angular speed of Y/frame. The histogram in Fig. 5B reproduces the average amount of the estimated expansion, rotation, and two components of shear. The histograms of Figs. 5C and D compare respectively the speed and direction distributions of the optical flow of Fig. 5A and of the corresponding noise free optical flow. As in the previous example, it is easy to realize that the algorithm is able to compute both optical flow and the elementary motions correctly and is almost not sensitive to the added noise. To conclude, let us mention a few differences of this method from previous methods. First, the present algorithm does not require the solution of differential equations, nor global functional minimization as in [24, 261. Second, optical flow is uniquely determined by two local assumptions: the changing image brightness is assumed to be stationary in time and optical flow can be approximated locally by a linear two-dimensional vector field.

Thus the computation of higher order derivatives of the image brightness pattern is not required (see [lo, 27-31]), because first-order derivatives are considered over a neighborhood of the point of interest. This fact means that an optical flow can be obtained which shows a high degree of spatial coherence and does not need further refining stages (see [lo], for example). Intuitively, the assumption of local linearity acts like a “shape regularization” constraint. It has to be noted, however, that from the qualitative point of view most of the various methods which have been proposed in the computer vision literature for the computation of optical flow produce similar results [32]. This is mainly due to the simple spatial structure of optical flow of rigid objects discussed in Section 2. 5.

3D MOTION

RECOVERY

In this section, the first order properties of optical flow are effectively used to distinguish between qualitatively

98

CAMPANI AND VERRI

80 60

a

60

60

0

0

P

Y

50 60 70 80

6

90

loo

FIG. 5. (A) The optical flow obtained through the proposed method while the noisy pattern of Fig. 3 was rotating around an axis orthogona1 to the image plane (N = 41 and n = 42). (R) Average estimates of the amount of expansion, rotation, and two components of shear. (C) Direction distribution and (D) Speed distribution for the optical flow obtained from the noise-free sequence (white columns) and the optical Bow of (A) (textured columns). The histograms are normalized to the maximum value.

different kind of motions and estimate 3D motion param- collision, and the polygonal line, which accounts for the eters. estimated time-to-collision, is typically below 10%. From Figure 6A shows one frame of a sequence in which the Fig. 6D it is evident that the time-to-collision has often viewing camera is translating toward a picture posted on been overestimated. This is mostly due to the average the wall. Figure 6B reproduces the optical flow computed underestimation of the optical flow components caused by means of the proposed method for the frame of Figure by the smoothing step which is required to compute im6A, while the average amount of the estimated expan- age brightness derivatives. Similar results have been obsion, rotation, and two component of shear is displayed tained on other sequences. To sum up, it appears that the in Fig. 6C. It is evident that both the flow field of Fig. 6B reconstructed optical flow is usually sufficiently accurate and the computed elementary motions of Fig. 6C agree to allow for the recognition of translational motion and qualitatively with the expanding flow field that is ex- the estimation of the time-to-collision. pected from a pure translation of the viewing camera Let us briefly describe the procedure which has been toward a Aat surface parallel to the image plane. The adopted for the computation of the time-to-collision. large percentage of expansion with respect to the other First, the focus of expansion p is located as the point elementary motions can be used to identify the viewed which minimizes the sum of the distances from the motion as “translational motion. ” straight lines that go through each flow estimate. Then, Figure 6D shows the estimated time-to-collision over the linear vector field VIcentered in p which best approxithe entire sequence. Note that the intersection of the mates the estimated optical flow is computed. The timecoordinate axes does not coincide with the origin. In fact, to-collision is finally calculated as the average of the inthe agreement between the straight line, the true time-to- verses of the diagonal elements of the matrix which

MOTION ANALYSIS

FROM OPTICAL FLOW

99

C

80 60 40 20

0 FIG. 6. (A) A frame of a sequence in which the camera (Pulnix TM46) was translating toward a picture posted on the wall. (B) The corresponding optical flow computed by means of the described method (N = 41 and n = 42). (C) Average estimates of the amount of expansion, rotation, and two components of shear for the optical flow of(B) (normalized to the maximum value). (D) Plot of the estimated (solid line) and true (dashed line) time-to-collision against time ([frames]) over the entire sequence.

defines the linear vector field vi. This procedure appears to be rather robust (given the large number of flow estimates, that is, the redundancy in the computation of the flow field) and able to locate the focus of expansion even outside the field of view (that is, when the direction of motion is not orthogonal to the image plane). Let us now discuss a second example. Figure 7A shows one frame of a sequence where a plant is moving on a rotating support. Figure 7B reproduces the estimated optical flow for the frame of Fig. 7A. The average amounts of the estimated expansion, rotation, and shears are displayed in Fig. 7C. From Fig. 7C it is again easy to identify the viewed motion as a “rotation” around an axis orthogonal to the image plane. The quantitative agreement of the reconstructed optical flow and elementary motions can be checked by computing the angular velocity. As can be easily seen in Fig. 7D, the estimated angular velocity, over the entire sequence, is very close to the true angular velocity. Let us now describe the procedure adopted for the

computation of the angular velocity (which is very similar to the procedure employed in the previous example). The center of rotation p is found as the point which minimizes the sum of the distances of the straight line that are orthogonal to each flow estimate through the base point. Standard least-mean-square techniques are then used to produce the linear vector field vI centered in p which best approximates the estimated optical flow. The angular velocity is finally estimated as the average of the difference of the antidiagonal elements of the matrix which defines the linear vector field VI. In conclusion, the proposed method proved that it can reliably distinguish between different motions and estimate motion parameters quite accurately. 6.

SURFACE

RECONSTRUCTION

In this section it is shown that the first-order properties of optical flow can be helpful in estimating surface depth and orientation in the case of camera translation.

CAMPANI AND VERRI

100

:::::

a

2

3

4

5

6

7

FIG. 7. (A) A frame of a sequence in which a plant was rotating around an axis nearly parallel to the camera optical axis. (B) The corresponding optical flow computed by means of the described method (N = 41 and n = 42). (C) Average estimates of the amount of expansion, rotation, and two components of shear for the optical flow of(B) (normalized to the maximum value). (D) Plot of the estimated (solid line) and true (dashed line) angular velocity ([“/frame]) against time ([frames]) over the entire sequence.

Figures 8A, B, C, and D show four frames of four sequences in which the viewing camera is translating toward the same, but differently oriented, planar surface (the angle cpbetween the unit normal to the viewed surface and the direction of motion equals o”, W, 30”, and 45”, respectively). The corresponding optical flows, computed through the method described in Section 4, are shown in Figs. 9A, B, C, and D. Let us now illustrate two simple methods which can be used to distinguish between the rather similar optical flows of Figs. 9 and recover surface orientation and depth. The first method is based on the simple formula which relates the position of the focus of expansion x0 and the optical flow v at each location x of the image plane with the perceived depth T(X) at that location; that is,

T(x)= 11x ,$ v

’

(11)

where G-(X)in Eq. (11) is actually the time-to-collision and can be used to estimate absolute depth if the viewer

translational velocity is known. Table 2 shows the reconstructed orientation y31and distance rl of the planar surfaces of Fig. 8 which are obtained by computing the equation of the planar surface that best approximates the depth estimates obtained from the optical flows of Fig. 9 through Eq. (11) at each location where optical flow is TABLE 2 Estimates of Orientation and Distance of a Translating Plane (0

7

X7

-d

PI

0 15 30 45

100 104 115 130

127 127 124 125

126 125 130 142

o-c3 14 ? 4 29 t 3 45 + 3

71 102 2 106 f 116-t136 f

w 1 1 1 1

6k6 22 t 4 29 k 5 50 2 6

72 58 c 100 + 121 * 140 +

8 4 5 9

Nofe. cpand r are the true values of the angle between the optical axis and the normal direction to the plane ([“I) and of the time-to-collision ([frames]) respectively, while xy and xi are the estimated components of the singular point for each image sequence, p,, 71 and ppz,72are the values computed through the first and the second method respectively described in the text. The uncertainties are the standard errors.

MOTION ANALYSIS

FROM OPTICAL FLOW

101

FIG. 8. (A), (B), (C), and (D) Four frames of four image sequences in which a differently oriented planar surface is translating toward the image plane (the angle p between the unit normal to the viewed surface and the direction of motion equals O”, 15”, 30”, and 45” respectively).

computed. The focus of expansion is located according to fTIT3 is the location of the focus of expansion on the the procedure described in Section 5. It is evident that image plane. From Eq. (12) it is easy to compute the this method is able to qualitatively distinguish between entries of the matrix M which read the different orientation of the viewed planar surfaces. Quantitative estimates are also very good. It is worthwhile mentioning that, due to the small aperture of the viewing device, in the experimental conditions in which 1 XI - x: the method has been tested, an error of 1% in the depth M12 = ; ~ U’X u2 estimate (nearly corresponding to 1 cm) is responsible for an error of about 2” in the reconstructed slant. (13) 1 x2 - x!j M2, = - ~ Ul Let us now describe the second method which is based 7 u-x on an analysis of the elementary motions of optical flow. The optical flow of a translating, arbitrarily oriented, planar surface which at time t satisfies the equation u * X = g(t) can be written as where 7 = -X3/T3 is the time-to-collision of the point X (that is, as before, the distance in time units of the point u . x(x - x”), (12) X on the planar surface from the image plane). Equations ” = - j-:t, (13) can be rewritten as where T3 is the component of the translational velocity bi.u=O i=l,. . . 74, T along the optical axis, f is the focal length, and x0 = (14)

CAMPANI AND VERRI

102

FIG. 9. (A), (B), (C), and(D) The optical flows associated with the sequences of Fig. 8 and computed through the method described in Section 4 (N = 41 and n = 17).

where Axi = xi - xp, i = 1, 2. In the specific example of Fig. 8, Eq. (16) can be simplified by imposing the constraint that the viewed surface is orthogonal to the optical axis. Since this implies that ~2 = 0, Ml2 = 0, and

where the vectors bi are bl = (TM~, - 1)x - (xl - x$e, b2 = rM12x - (x1 - xy)e;! b3 = 7MZIx - (x2 - xi)e,

(15)

b4 = (7Mz2 - 1)x - (x2 - xi)e2.

That bl x b3 = 0 and b2 x b4 = 0 can be easily checked. Therefore, the bi in Eq. (15) identify a plane in 3D space. From Eqs. (14) it can be concluded that this plane is parallel to the viewed planar surface. The direction parallel to the unit vector u can thus be estimated by taking the vector product of b2 and b3. Explicit computation yields b2 x b3 = T(Ax1M2,e2 - Ax2M12el) x x + AxlAx2e3, (16)

M22 = l/7,

(17)

Eq. (16) can be rewritten as b2 x b3 = Ax,

2

(18)

Equation (18) has been used to estimate the orientation of the planar surface at each point where optical flow is computed. The results which have been obtained by applying this second method to the image sequences of Fig. (8) are also shown in Table 2. The estimated orientation,

MOTION

ANALYSIS

‘p2, has been determined by averaging the orientation obtained at each patch in which optical flow has been computed. The surface distance is measured as the time-tocollision, which can be obtained from Eq. (17), where the matrix M is computed at the patch that contains the focus of expansion. In order to improve accuracy, the method requires the computation of optical flow over larger patches (in this case N = 81 and IZ = 9). It is evident that the method is sufficiently accurate to make it possible to qualitatively distinguish between the differently oriented surfaces. Quantitatively, however, this method appears to be less reliable than the previous technique. This is probably due to the fact that the information which is used in this second case is entirely contained in the differential property of optical flow, as opposed to the first method in which only optical flow estimates are used. 7. IMAGE SEGMENTATION

In both of Sections 5 and 6, the presented optical flow analysis relies upon the assumption that the flow estimates refer to the same surface. This assumption is usually unrealistic in the absence of a method which is able to segment the viewed image into different surfaces. In this section an extension of a method recently proposed for image segmentation [33] is described. The method, which makes use of the first-order spatial properties of optical flow, consists of three main steps. In the first step the image plane is divided into a number of nonoverlapping squared patches (typically of 16 x 16 pixels in 256 x 256 images) and first-order properties of optical flow are analysed. In the second step, a list of the “possible” viewed motions, or labels, is produced. In the third step, according to a suitable similarity measure, a label, that is, a possible motion, is attached to each patch by means of a relaxation technique. It should be noted that the idea of labeling image patches depending on the apparent motion was first proposed in [34, 351. A thorough discussion of the main similarities and differences between the presented method and the method described in [35] can be found elsewhere [36]. Figure 10A shows a frame of a synthetically generated sequence in which the larger sphere is translating toward the image plane, while the smaller sphere is translating outward the image plane and the background is rotating. The optical flow relative to the frame of Fig. lOA, and computed through the procedure described in Section 4, is shown in Fig. 10B. In order to identify the different motions in Fig. lOB, the first step of the method analyzes the first-order spatial properties of optical flow. The optical flow is divided into patches of fixed size and the expanding (EVF), contracting (CVF), clockwise (CRVF) and anticlockwise (ARVF) rotating, and constant (TVF) vector fields which best approximate the optical flow in

FROM OPTICAL

FLOW

103

each patch si are computed. Roughly speaking, this is equivalent to reducing the possible 3D motions to translation in space with a fairly strong component along the optical axis (EVF and CVF), rotation around an axis nearly orthogonal to the image plane (CRVF and ARVF), and translation nearly parallel to the image plane (TVF). This choice, which is somewhat arbitrary and incomplete, does not allow for an accurate recovery of 3D motion and structure (the shear terms, for example, are not taken into account), but usually seems to be sufficient to obtain a qualitative segmentation of the viewed image in the different moving objects [36]. As a result of the first step, five vectors x< , j = 1, . . . , 5, are associated with each patch si: the”vector xi, position over the image plane of the focus of expansion of the EFV, x2, position of the focus of contraction of the CVF, xi,, pos%ion of the center of the CRVF, x:~, position of the center of the ARVF, and the unit vector x5s,’ parallel to the direction of the TVF. In order to produce a list of the “possible” motions, in the second step, global properties of the obtained EVFs, CVFs, CRVFs, ARVFs, and TVFs are analyzed. This step is extremely crucial, since the pointwise agreement between each of the computed local vector fields and the optical flow of each patch usually makes it impossible to select the most appropriate label. Figures 1OC and D, respectively, show the distribution of the foci of expansion and contraction and centers of clockwise and anticlockwise rotation associated with the EVFs, CVFs, CRVFs, and ARVFs of the optical flow of Fig. 10B. A simple clustering algorithm has been able to find two clusters in the distribution of Fig. IOC, and these clusters clearly correspond to the expansion and contraction along the optical axis of Fig. 10B. The same algorithm, applied to the distribution of the centers of rotation (shown in Fig. lOD>, reveals the presence of a single cluster in the vicinity of the image plane center corresponding to the anticlockwise rotation in Fig. 10B. On the other hand, in the case of translation, the distribution of the unit vectors parallel to the directions of the TVFs is considered (see Fig. IOE). For the optical flow of Fig. 10B the distribution is nearly flat (see Fig. 10E) indicating the absence of preferred translational directions. Therefore, as a result of this second step, a label 1is attached to each “possible” motion which can be characterised by a certain cluster of points x:“, where c(l) equals 1, . . . , 4, or 5 depending on 1. In the specific example of Fig. 10, one label of expansion, one of contraction, and one of anticlockwise rotation are found. In the third and final step, each patch of the image plane is assigned one of the possible labels by an iterative relaxation procedure (deterministic version) [37]. The key idea is to define a suitable energy function which depends on both the optical flow patches and the possible

104

CAMPANl AND VERRI

c

D

384

384

320-

. +

r.’

1

320

1

* l

256 .

192 1

128

*

*

l +

+++

i

:

l

+.

*+i*

t

*

#+‘*‘ .

.

0

2

+

+ *+

64

1

*

t+ *

+

,

t *,*+ .

**

l

I

;

-64 -128

-64

0

64

128

192

256

320

384

-128

-> 12 8-64

l

0

1

64

128

192

256

320

384

F

FIG. 10. (A) A frame of a synthetic sequence in which the larger sphere is translating toward the image plane, while the smaller sphere is moving away and the background is rotating anticlockwise. (B) The corresponding optical flow computed by the method proposed in Section 4. (C) Distributions of the foci of expansion (squares) and contraction (crosses) of the EVFs and CVFs, respectively, which lie within an area four times larger than the field of view (identified by the solid frame). (D) Distribution of the centers of anticlockwise rotatio of the ARVFs. (E) Distribution of the directions of the TVFs on the unit circle. (F) Color coded segmentation of the optical flow of(B) obtained through the algorithm described in the text.

motions and reaches its minimum when the correct Iabels are attached to the flow patches. In the current implementation, the energy function is a sum extended over each pair of neighboring patches in which the generic

term u(s~, si), where S; and sj are a pair of neighboring patches, is given by the formula

MOTION ANALYSIS

FROM OPTICAL FLOW

105

C 384, 320. 256. 192. 128. 64 O-64. -1281 -128

’

4

-64

0

64

128

192

256

320

384

FIG. 11. (A) A frame of a sequence in which the box was translating toward the camera, while the camera was translating toward an otherwise static environment. (B) The corresponding optical flow computed by the method of Section 4. (C) Distribution of the foci of expansion of the EVFs. (D) Color coded segmentation of the optical flow of (B) obtained through the algorithm described in the text.

where xI is the center of mass of the cluster corresponding to the label 1, and 6 = I if the two patches have the same label, otherwise 6 = 0. The relaxation procedure has been implemented through an iterative deterministic algorithm in which, at each iteration, each patch is visited and assigned the label which minimizes the current value of the energy function, keeping all the other labels fixed. The procedure applied to the optical flow of Fig. lOB, starting from a random configuration, produces the color coded segmentation shown in Fig. IOF after 20 iterations. From Fig. lOF, it is evident that the method is able to detect and correctly identify the multiple motions of the optical flow of Fig. 10B. Let us finally discuss an example of real images. Figure 11A shows a frame of a sequence in which the viewing camera is translating toward the scene while the box is moving toward the camera. The optical flow associated with the frame of Fig. 11A is shown in Fig. 11B. From Fig. 11B it is evident that the problem of finding different moving objects from the reconstructed optical flow is not

simple. Due to the large errors in the estimation of optical flow, simple deterministic (and local) procedures which detect flow edges, or sharp changes in optical flow, are doomed to failure. In addition, the viewed motion consists of two independent expansions and even in the presence of an exactly computed optical flow, no clear flow edge can be found in the vicinity of the top, right side, and bottom of the box. Figure 11C shows the distribution of the foci of expansion associated with the EVFs computed as described above. Two clusters are found which correspond to the (independent) motion of the camera and of the box of Fig. 11A. In contrast, no clusters are found in the other distributions. Therefore, it can be concluded that, at most, two different motions (mainly along the optical axis) are present in the viewed scene. The color-coded segmentation which is obtained by applying the third step of the proposed method is shown in Fig. 11D. It is evident that the algorithm detects and correctly identifies the two different motions of the viewed scene. To summarize, the present method seems to be able to

106

CAMPANI

detect multiple motion and correctly segment the viewed image in the different moving objects even from rather noisy and imprecise estimates of optical flow. Due to the small number of labels, or possible motions, and to the coarse resolution at which the segmentation step is performed, the algorithm only takes a few seconds on a Sun SPARCStation on a 256 x 256 image, apart from the computation of optical flow.

AND VERRI

quence of images, in Proceedings,

3rd Internat. Vision, Osaka, Japan, 1990, 441-445.

6. J. J. Koenderink and A. J. Van Doorn, How an ambulant observer

8. 8. CONCLUSIONS

This paper has shown that the first order properties of optical flow, such as the singular points and elementary motions, are usually sufficient to describe and analyze image motion. A number of conclusions can be drawn from our analysis. 1. The optical flow of rigid objects, apart from surface and motion boundaries, can usually be approximated, over rather large areas of the image plane, by a linear vector field. 2. Optical flow based methods can rely upon patchwise, instead of pointwise, flow estimates. This makes it possible to devise algorithms for the analysis of moving images which are less sensitive to noise. 3. The singular points and the elementary motions can be used to distinguish between the different kind of motions, estimate motion parameters, extract information on surface depth and orientation, and segment optical flow into the different moving objects. The intrinsic spatial coherence of patchwise estimates of optical flow makes it possible to solve these problems effectively. ACKNOWLEDGMENTS This work has been partially funded by the Progetto Finalizzato Robotica, the PROMETHEUS project (Progetto Finalizzato Trasporti), by the ESPRIT project VOILA, and by the Agenzia Spaziale Italiana. M. C. has been partially supported by the Consorzio Genova Ricerche. Clive Prestt checked the English.

REFERENCES

Conf. Comput.

9. 10.

11.

12.

13.

can construct a model of the environment from the geometrical structure of the visual inflow, in Kibernetic 1977 (G. Hauske and E. Butendant, Eds.), Oldenbourg, Munchen, 1977. K. Prazdny, Determining the instantaneous direction of motion from optical flow generated by a curvilinear moving observer. Comput. Vision Graphics Image Process. 22, 1981, 238-248. L. Dreschler and H. H. Nagel, Volumetric model and 3D trajectory of a moving car derived from monocular TV frame sequences of a street scene, in Comput. Vision Graphics Image Process. 20, 1982, 199-228. J. J. Koenderink, Optic flow, Vision Res. 26, 1986, 161-168. S. Uras, A. Verri, F. Girosi, and V. Torre, A computational approach to motion perception, Biol. Cybernet. 60, 1988, 79-87. A. Verri, F. Girosi, and V. Torre, Mathematical properties of the two-dimensional motion field: From singular points to motion parameters, J. Opt. Sot. Am. A 6, 1989, 698-712. R. C. Nelson, and J. Aloimonos, Obstacle avoidance using flow field divergence, IEEE Trans. Pattern Anal. Mach. Intelligence 11, 1989, 1102-l 105. G. Sandini, M. Tistarelli, Robust obstacle detection using optical Bow, in Proceedings, First IEEE Workshop Robust Comput. Vision, Seattle,

WA, 1990.

14. J. J. Gibson, The perception of the visual world, Houghton Mifflin, Boston, 1950. 15. R. Y. Tsai, T. S. Huang, and W. L. Zhu, Estimating 3D motion parameters of a rigid planar patch. II. Singular value decomposition, IEEE Trans. Acoust. Speech Signal Process. 30, 1982. 16. H. C. Longuet-Higgins, The visual ambiguity of a moving plane, Proc. R. Sot. London R 223, 1984, 165-175. 17. S. J. Maybank, The angular velocity associated with the optical flow field due to a single moving rigid plane, in Proceedings, 6th European Conf. Arttjicial Intell., 1984, pp. 641-644. 18. S. Negahdaripour and B. K. P. Horn, Direct passive navigation, IEEE Trans. Pattern Anal. Mach. Intelligence 9, 1987, 168-176. 19. M. W. Hirsh and S. Smale, Differential Equations, Dynamical Systems, and Linear Algebra, Academic Press, New York, 1974. 20. A. Malisia, M. Campani, and A. Verri, A simple and robust method for estimating the motion of a mobile vehicle from optical flow, submitted for publication. 21. H. Helmholtz, Uber Integrale der hydrodynamischen Gleichungen welche den Wirbelwegungen entsprechen, Crelles J. 55, 1858, 25. 22. M. Campani and A. Verri, Computing optical flow from an overconstrained system of linear algebraic equations, in Proceedings, 3rd Internat. Conf. Comput. Vision, Osaka, Japan, 1990, pp. 2226.

1. R. Y. Tsai and T. S. Huang, Uniqueness and estimation of 3D motion parameters of rigid objects with curved surfaces, IEEE Trans. Pattern Anal. Mach. Intelligence 6, 1984, 13-27. 2. A. M. Waxman, Image flow theory: A framework for 3-D inference from time-varying imagery, in Advances in Computer Vision (C. Brown, Ed.), Erlbaum, New Jersey, 1986. 3. 0. Faugeras, F. Lustman, and G. Toscani, Motion and structure from point and line matches, in Proceedings, 1st Internat. Conf. Comput. Vision, London, U.K., 1987, pp. 25-34. 4. T. Vieville and 0. Faugeras, Feed-forward recovery of motion and structure from a sequence of 2D-lines matches, in Proceedings, 3rd Internat. Conf. Comput. Vision, Osaka, Japan, 1990, pp. 517-520. 5. J. L. Jezouin and N. Ayache, 3D structure from a monocular se-

23. C. L. Fenmena and W. B. Thompson, Velocity determination in scenes containing several moving objects, Comput. Graphics Image Process. 9, 1979, 301-315.

24. B. K. P. Horn and B. G. Schunck, Determining optical flow, Artif. Intelligence

17, 1981, 185-203.

25. A. Verri and T. Poggio, Motion field and optical flow: Qualitative properties, IEEE Trans. Pattern Anal. Mach. Intelligence 11, 1989 490-498. 26. E. C. Hildreth, The Measurement

of Visual Motion, MIT Press, Cambridge, MA, 1984. 27. H. H. Nagel, Displacement vectors derived from 2nd order intensity variations in image sequences, Comput Vision Graphics Image Process. 21, 1983, 85-I 17.

MOTION ANALYSIS 28. R. M. Haralic and J. S. Lee, The facet approach to optic flow, in Proceedings, Image Understanding Workshop (L. S. Baumann, Ed.), pp. 84-93, Science Applications, Arlington, VA, 1983. 29. 0. Tretiak and L. Pastor, Velocity estimation from image sequences with second order differential operators, in Proceedings, International Conference on Pattern Recognition, Montreal, 1984, pp. 16-19. 30. E. H. Adelson and J. R. Bergen, Spatiotemporal energy models for the perception of motion, /. Opt. Sot. Am. A 2, 1985, 284-299. 31. A. Verri, F. Girosi, and V. Torre, Differential techniques for optical flow, J. Opt. Sot. Am. A 7, 1990, 912-922. 32. J. J. Little and A. Verri, Analysis of differential and matching methods for optical flow, in Proceedings, IEEE Workshop on Visual Motion, Irvine, CA, 1989, pp. 173-180.

FROM OPTICAL FLOW

107

33. A. Rognone, M. Campani, and A. Verri, Detecting moving objects from optical flow, Pattern Recognit. Image Anal. 2(l), 1992. 34. P. Bouthemy and J. Santillana Rivero, A hierarchical likelihood approach for region segmentation according to motion-based criteria, in Proceedings, 1st Internat. Conf. Comput. Vision, London, 1987, pp. 463-467. 35. E. Francois and P. Bouthemy, Derivation of qualitative information in motion analysis, Image Vision Comput. 8, 1990, 279-288. 36. A. Rognone, M. Campani, and A. Verri, ldentifying multiple motions from optical flow, in Proceedings, 2nd European Conference on Computer Vision, St. Morphenite Lipure, Italy, 1992, pp. 258266. 31. D. Geman, S. Geman, and C. Graffigne, and P. Dong, Boundary detection by constrained optimization, IEEE Trans. Pattern Anal. Mach. Intelligence 12, 1990, 609-628.