- Email: [email protected]
Structure-from-motion algorithms for computer vision on an SIMD architecture
Computer Physics Communications 37 (1985) 273—280 North-Holland, Amsterdam
273
STRUCTURE-FROM-MOTION ALGORITHMS FOR COMPUTER VISION ON AN SIMD ARCHITECTURE B.F. BUXTON and D.W. MURRAY GEC Research Laboratories, Hirst Research Centre, East Lane, Wembley, UK
H. BUXTON and N.S. WILLIAMS Department of Computer Science and Statistics, Queen Mary College, Mile End Road, London El 4NS, UK
Algorithms for interpreting the motion of edge features in an image sequence in terms of the position, orientation and motion of a planar visible surface facet have been developed and implemented on an SIMD processor array. The underlying theory of the interpretation based on a simultaneous solution of the structure-from-motion and aperture problems is briefly reviewed and it is explained how the solution can be implemented very simply and efficiently on an SIMD processor array. Test calculations were carried out on the ICL DAP at Queen Mary College to illustrate the implementation and to show that, for an ideal 1 : 1 mapping of the image pixels onto the processor array, the algorithms run at speeds that are nearly compatible with real time video frame rates.
1. Introduction: computer vision Most of us are familiar with the idea that depth perception can be attained by stereo vision. Each eye (or sensor) sees a different view of the world with the result that the images formed are slightly out of alignment. From this shift or ‘disparity’ the depth can be computed with the help of a little geometry. However, it is not so widely appreciated that depth information can also be obtained if the sensor and the objects in view are moving relative to each other. One way to see this is to regard image snapshots taken at successive positions as stereo pairs. The two problems, both involving reconstruction of the depth dimension ‘missing’ in an ordinary image are obviously related. Moreover, in both cases it is the interpretation of the images in .terms of the three dimensional scene in view that generally distinguishes these problems as computer vision rather than belonging to the realm of image processing and pattern recognition. In image processing or pattern recognition the image is analysed as a two dimensional pattern per se, even when it is the projection of three dimensional objects in a three dimensional world [1].
The application of vector and parallel processors to image processing is now fairly commonplace, especially if SIMD or pipeline machines are used, and a number of more complicated architectures have been suggested for pattern recognition tasks (see ref. [21 for a recent review). Application of parallel computers in computer vision is not so widespread, but SIMD machines have been used for some low-level tasks and MIMD structures have been proposed for higher level vision work. Brady [3]has recently given a survey of parallelism in computer vision, including reference to some dedicated hardware implementations. In this paper, we want to describe the algorithms we have recently derived for obtaining 3D information from visual motion [4] and discuss their implementation on one of the simplest parallel architectures, an SIMD machine. The results of test calculations on the ICL DAP at Queen Mary College [5—7]will be used to illustrate the details of the implementation. The paper is organized as follows. We begin in the next section with a brief description of how optic flow (the technical term used to describe the motion of points in the image plane of the sensor)
OO1O-4655/85/$03.30 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
B.F Buxton et a!. / Structure-from - motion algorithms on an SIMD architecture
274
contains information about the depth and orientation of the visible surfaces and their motion relative to the sensor. We also explain why it is convenient to compute the motion of edge features in an image sequence rather than the whole of the optic flow field. As this only gives partial information about the flow field, we have to solve an aperture problem arising from the lack of information about movement of the edge features along their length in addition to reconstructing the missing depth information. The solution obtained recently [4] is outlined in section 3 and implementation of the algorithms derived is described in section 4. Particular attention is paid to the very simple way in which these algorithms map onto an SIMD processor array and the very high processing speeds that can be achieved. Geometrical conditions that cause the algorithms to fail and the numerical stability of the calculations are also briefly discussed. We end with a short conclusion and discussion that includes some comments on the direction of future work.
motion over the whole of the image plane have been discussed in a number of papers [8—12],and it has been shown that, in principle, the optic flow field contains enough information to allow one to use it to reconstruct the structure of rigid bodies moving relative to the sensor [13,14]. Obtaining the structure of a rigid body this way is referred to as the structure-from-motion problem. In our work, we consider a slight variant of the structure-from-motion problem in which we attempt to reconstruct visible surfaces from the visual motion. To show that the optic flow field contains information about the position and orientation of the visible surfaces, it is convenient to assume that the object points, R, lie on a planar surface facet, R N 1 2 3 — —
moving (instantaneously) with rectilinear velocity V and angular velocity ~Q with respect to the sensor. The optic flow field, r(r), given by 1=
V(rN)+r(~ V)(rN)/l+Q/\r +r[~.(QAr)]/l,
(2.4)
2. Optic flow As mentioned above, optic flow is the motion of points in the image plane due to motion of the corresponding object points relative to the sensor. The velocity of an image point r depends on the position, R (X, Y, Z), and velocity, R, of the object point according to (2.1) t 1R/Z+ IRZ/Z2, =
if we use the perspective transformation =
—
IR/Z
(2.2)
to model the geometry of the sensor (fig. 1). Details of the optic flow field, r(r), describing the
I
~J Fig. 1. The sensor geometry.
________>~
then depends on the surface and motion parameters N, V and £2. Since the terms involving N and V in (2.4) always appear as products of their components, we note that there is an inherent ambiguity in the optic flow field as it is the same for a small object close to, moving slowly, as for a large object, far away, moving quickly. Given the optic flowber at a small number of image points, (2.4) can solved for the surface and motion parameters [15—17]to within this inherent ambiguity (and a few other complications that we shall mention later in section 3). Thus one c,,an determine the direction of rectilinear motion V. the surface orientation, N, the angular velocity £2 and the the value of the products VN or v-N [15]. However, before these algorithms can be used, observed dE~)O th~ doing thisimage is to follow the motion of prominent image features from frame-to-frame, but this involves identifying prominent features in one frame and finding the corresponding features that they
B.F. Buxton et at
/
Structure-from - motion algorithms on an SIMD architecture
match in later frames. Algorithms of this kind involving a solution to the correspondence problem [18,19] are closely related to binocular stereo algorithm [20—24].There is however a simpler alternative approach, based on the motion constraint equa(ion [25—26], 3E ‘at —r E 2 ~ —
/
—
~‘
‘
/
‘
that avoids the correspondence problem. Since stands for the gradient operator (a/ax, a/ay) in the image plane, eq. (2.5) implies that the total derivative d E/d t is zero. Using the motion constraint equation is therefore equivalent to following the motion of a given image feature point which is usually assumed to be the same as obtaining the optic flow defined geometrically above, A fundamental difficulty with (2.5) however is that it can only yield the component of the optic flow parallel to ~E denoted here by r~as it is perpendicular to the edge contour. This lack of information about motion perpendicular to ç7E is known as the aperture problem, but the way it manifests itself depends on how (2.5) is used. Since the image irradiance values E ( r, t) are always corrupted by sensor noise, smoothing or filtering is required before (2.5) can be used. In addition, (2.5) will only give reliable estimates of i~where both aE/at and lyE I are large as is the case at moving edge features in the image [26—28].As described elsewhere [28,29], we therefore smooth the image irradiance values E(r, I) over both space and time with a Gaussian filter G(r, t) ~ 2~2±~2t2),~’2a2 and apply (2.5) at the edge feae_ tures defined by the zero-crossings of the visual signal [28,29], ~‘
S(r,
t)*E(r,
t)= _(~2+i~)G(r,
275
out the spatio-temporal convolution in (2.6), detecting the zero-crossings in S and estimating the vernier velocities v from (2.7), have already been described [28—30].We need only note here that as similar operations are to be carried out at every pixel in the image, these calculations map very simply onto an SIMD processor array such as the DAP. In addition, the masking facilities available on the DAP enable us to write simple uniform algorithms that only calculate the vernier velocity v at the edge features where it is required [30]. Linear least-squares interpolation enables us to calculate both v and the positions of the edge features to pub-pixel accuracy. As shown in fig. 2 for a toy calculation simulating a bright circular blob painted on a plane moving directly towards the sensor, the vernier velocities then accurately depict the edges of features on the visible surface. Other examples are given in the refs. [4,28—30]. The final point to note here is that these low-level calculations are extremely fast on the DAP, even when not fully optimized [28,30]. For the ideal 1: 1 mapping of the image pixels directly into the processor array, calculations of the kind illustrated in fig. 2 can be carried out in approximately 35 ms per frame, even though with the masks of radius 8 used (o 2%/L u I in the Gaussians), we are =
=
\\\
~ 1111 /f/ / / 1/,
.-
:
t).
(2.6) This yields an average value of the optic flow normal to the edge features, called the vernier velocity v [28]:
~
/
\“
~1, ~=~=
[(_as/at)(vs)] V
.
(2.7)
11,1 !•l~ ~ . Fig. 2. The vermer velocities at the edge of a brigbt circular
carrying
blobsensor. painted on the surface of a plane moving directly towards the
edge
Details of this calculation which involves
276
B. F. Buxton et a!.
/
Structure -from - motion algorithms on an SIMD architecture
effectively performing convolutions over 2000 pixels [30,31].
3. Interpreting edge motion The problem then is to interpret the edge motion by simultaneously solving the aperture and structure-from-motion problems. Eq. (2.4) relates the optic flow r to the motion and surface parameters V, £2 and N. To eliminate the unknown cornponent of the optic flow parallel to the image edge features, we take the dot product of (2.4) with the vernier v, leading to theV)(r.N)/l equation [4], 2 (v.velocity V)(r.N) +(vr)(z1 v +v.(S2Ar)+(v•r)[~(~2Ar)]/1, (3.1) =
for V, £2 and N. By writing the vectors in terms of their Cartesian components, V (V 1, V2, V3) etc., and using the summation convention, (3.1) may be rearranged in the form [4],
may be inverted for the parameters p,~leaving us to solve the eight non-linear simultaneous eqs. (3.4). The latter, however, is a problem that was encountered a little while ago by Longuet-Higgins in his work on the interpretation of the image of a moving planar surface [15], and we can take over the results of his analysis. Briefly, the solution proceeds by eliminating the unknowns that appear linearly in (3.4), the ~k, and constructing a cubic equation for the ‘odd-man-out’ amongst the remaining quantities, 173 N3. Once 173 N3 has been obtained from this cubic, we can solve for the remaining V ]VJ. and finally for the ~ Note, as mentioned 2, that because only 1’J~occur ininsection (3.4), the magnitudes of Vproducts and N cannot be fixed independently. In addition, the vectors V and N can be interchanged and their signs altered without affecting the value of the products 1’~N~giving four interpretations of the
=
[v1~~
=
—(v.r)~e3x1/x31p,1.
(3.2)
where we have introduced the parameters P,i (3.3) and defined x3 to be —l so that r (x1, x2, x3). As it stands, (3.2) looks like an equation for the nine parameters p~,but since the vernier velocity v always lies in the image plane, v3 is zero with the result that only the combinations ~ p33 and P22 p33 appear in (3.2). We therefore rewrite (3.2) as a set of linear simultaneous equations for the eight parameters, —
=
—
—
Pj1
=
~
—
~ijk~k
—
V3N36~1,
(3.4)
at each of the edge feature points, m, in matrix form [4]: 2(m). (3.5) ~ c~,(m)p11 v =
visual motion. Although only one of these would be consistent with the observed edge feature positions at later times, information immediately at hand can be used to reduce the number of possible interpretations. This is because we know that the surface facet described by N must be visible from each of the feature points R(m) used in the solution. R(m).N must therefore be positive leaving, at most, only two viable interpretations. If in addition, there are points m such that R(m) V> 0 and other points m’ for which R (rn’)’ V < 0, only one interpretation is consistent with all the data. Thus, Longuet-Higgins [15] states the theorem that there is a unique interpretation if there are visible points in both the forward and backward fields of view with respect to the direction of rectilinear motion V. As Longuet-Higgins [15] has shown that the middle root of the cubic equation gives the required value of V3N3 (the others, if distinct, give rise to complex flow fields) the solution looks straightforward, identity for cos 30 especially can be since used to thecalculate trigonometric V
3N3 directly. The only remaining problem is the solution of the simultaneous eqs. (3.5) for the parameters p,~,which still involves processing the data at
The matrix c~(m)is defined by: (3.6) With data at eight or more edge points, m, (3.6)
the edge points m. In principle, this too is simple enough: we either solve the set of eqs. (3.5), which may be written in matrix form for n edge features,
B.F Buxton et a!.
/ Structure -from
-
motion algorithms on an SJMD architecture
as
(ii) calculate the components of the coefficient 2
(3.7) directly in a least squares sense, or use the classical least squares solution, [c]5~g[p]8xj[v
I
277
Ti
]nxS,
r i.....r T21
v~cJgxsl.pJ
—
1C
V
J8X1~
.
/
Since the equations are not usually very well conditioned, the former would be preferable, but we have not found a convenient routine available to calculate the pseudo-inverse of c on the DAP (e.g. by Householder transformations as in the NAG routines) when the matrix is stored with its long dimension (n) mapped onto a DAP plane and its short dimension (8) mapped vertically. We have therefore used the classical least squares solution (3.8) as the data is arranged almost ideally for this calculation (section 4).
4. Implementation on the DAP Before describing the details of the implementalion of our algorithms on the DAP at Queen Mary College, it is convenient to summarize the computational steps involved in calculating the motion and surface parameters, V, £2, N and (NV). For simplicity we shall assume at this stage that all visible points lie on the same planar facet. In general, of course, this will not be true and we would like to use optic flow to segment the scene into objects moving against the background, visible surface facets in various orientations, etc. [28,29]. Current work [32] suggests that an iterative segmentation scheme based on the simple calculation we shall describe below can be devised. We comment further on this in the discussion (section 5), but note here that the present calculation, in which it is assumed that there is a single set of surface and motion parameters valid across the whole image, may be of direct relevance to passive navigation [33]. When only a single set of surface and motion parameters are required, we have to: (i) obtain the position of the edge features r( m) and calculate the vernier velocities v( m) describing the edge motion,
2], matrix c defined (iii) compute 1(m) [cTc]and [cTvby (3.6), (iv) solve the least squares eqs. (3.8) for the parameters p,~, (v) find the middle root, V 3 N3, of the cubic
equation, (vi) use V3N3 to compute the surface and motion parameters N, V, £2 and (VN), (vii) use the position of the edge features to resolve the interpretations of V and N as far as possible. Note that steps (i)—(iii) and (vii) involve calculations over the image space (x, y) whereas (iv), (v) and (vi) involve only the small space of the eight surface and motion parameters. In the former therefore, it is most convenient to map the image pixels into the processor array. This has already been described [30] for step (i) involving convolution of the images as in (2.6), detection of the zero-crossing edge features and computation of the edge positions r( m) and vernier velocities v( m). A sequence of small 64 X 64 images can be mapped directly onto the DAP, but for larger images we have to map groups of several neighbouring pixels onto the same processor in a ‘crinkled’ or ‘pyramidal’ mapping scheme [30]. In either case however, the vectors describing the details of the edge position and the vernier velocities end up stored vertically below the processors in the DAP mem-
64 x 64 processor array
16K of memory
Fig. 3. Sketch of the architecture of the 64 x 64 ICL DAP showing the connectivity of the processing elements each haying access to its own memory.
278
B.F Buxton et a!.
/
Structure-from - motion algorithms on an SIMD architecture
ory planes (fig. 3). For simplicity, therefore, in what follows we shall assume the ideal 1: 1 direct mapping of the image pixels onto the DAP. Only small changes are required for larger images. We should note however that the masking facilities on the DAP are used to isolate the edge pixels so that no explicit encoding of the edge positions is required [30]. The r(m) are used to record the position of the zero-crossings to subpixel accuracy [28—30]. With the positions r( m) and the vernier velocities v( m) stored vertically below each edge pixel, it is then very easy to calculate the components of the coefficient matrix c,~(m)from (3.6) and store these vertically below each processor as well (fig. 4), where they are ideally positioned for computing [cTc] and [cTv2] in step (iii). This is because I cTc i ~ (~m)\ c, / \ ~ 1 ~m, — —
m
and 2] [cTv
(4.1) cj~(m) v2 (rn),
=
ni
enabling us to carry out the multiplications in parallel using planar multiplication on the DAP, and to perform the sums, ~m’ over the edge pixels using the planar sum function. The result is a very efficient routine that can be programmed in four or five lines of DAP Fortran [6,7], made even faster by the fact that cTc is, of course, a symmetric matrix. A statistical weighting w(m) of each point to improve the least squares fit as suggested by Scott [4,34,35] can also be incorporated with only minor changes.
_.-‘I~~’~ Columns of c
Processor array rows,m,of c”
,_‘7]
The next three steps, (iv)—(vi) in which the surface and motion parameters are computed could be carried out on an attached processor (the host machine in the case of the DAP), although in our test computations the least squares eqs. (3.8), step (iv), were also solved on the DAP in order to exploit the parallelism in the Gauss—Jordan elimination [6]. Steps (v) and (vi) were carried out on the host, however, and the results transferred back to the DAP for step (vii). Approximate times per frame for the steps carried out on the DAP are summarized in table 1, assuming as mentioned above the ideal 1: 1 mapping of the image pixels onto the processor array. All calculations were carried out in REAL *4 arithmetic except for the spatio-temporal convolution involved in step (i) which is done in 32-bit integer arithmetic (and takes roughly 20 ms, only half the interframe time at normal video rates). Note that a little time is ‘wasted’ on converting data from the DAP storage mode to the host 2980 storagebemode and back out again, operations which would better carried in separate hardware in a dedicated vision system [30]. Timings for steps (v) and (vi) on the host ICL 2980 are not given because in a vision system these calculations could be carried out on an attached serial processor (say) whilst the heavily pressed SIMD processor array continues with other calculations on steps (i) and (iii) (say). This kind of
Table 1 Approximate times per frame for the steps carried out on the DAP Step Function Time Machine (ms) (i) (ii)
(iv) (iii)
~ k L
(v) (vi) Fig. 4. Storage of the matrix c, 1(m) in DAP memory showing how row indices, m, map into the processor array and how the columns q and k! (say) are stored vertically in the memory planes.
(vii)
calculate r(m) and v(m) set up c. Tcl and [cTv2l 3(m) solve for the p,, calculate[c DAP-2980 storage mode conversion solve the cubic compute surface and motion parameters 2980-DAP storage mode conversion find correct interpretation
36
DAP
354 8 7
DAP DAP DAP
HOST HOST 3 4
DAP
B.F Buxton et a!. / Structure-from - motion algorithms on an SIMD architecture
concurrency would be possible, of course, since in applications video data would be continuously arriving for new processing to begin. Similarly, although we solved the least square equations (step (iii)), on the DAP, this could also be hived off to an attached processor as it only involves the small motion and surface parameter space. However, this question requires more careful consideration for a complicated scene when there would (after segmentation) be several sets of motion and surface parameters to be considered. Finally, we should consider the numerical stability of our algorithms. Solving for the surface and motion parameters from the p,~appears to be quite stable unless two roots of the cUbic for V~ N3 are coincident or nearly so. According to Longuet-Higgins’ analysis [15] this happens when V and N are parallel or antiparallel and the sensor is moving directly towards or away from the plane. Given the pt,, steps (v) and (vi) are thus usually quite stable but we remarked above and in section 3 that the least squares fit for the p~ is not well-conditioned. There are two reasons for this. The first is that all 3D scene interpretations derived from structure-from-motion algorithms seem to depend to some extent on perspective effects that will be small unless the points in view encompass a large visual angle [36]. This can be seen directly for our algorithm by referringback to eq. (3.1) and assuming that the image co-ordinates x and y are small compared to the camera length 1. Since r (x, y, —1), (3.1) reduces in lowest order to =
= =
IN3 ( v1V1 + v2V2) lv1Q2 + lv2~21 —1v1(V1N3 + £22) —1v2( V2N3 ~ —
—
—
(4.2)
showing that only the parameters P13 and P23 can be reliably determined if the field of view is small, This shows that all but two of the singular values of the matrix c.~(m)will be small as may be verified directly. It is also in agreement with our previous results [4] on the stability of these structure-from-motion algorithms when the motion parameters are known a priori and we have to solve (3.1) for the surface normal N. To see this we note that if V and £2 are known, we can use the values of
p13=V1N3+Q2
and
p23=V2N3—f2~
(4.3)
279
to estimate N3 accurately, as found in ref. [4]. The second reason why (iv) can be unstable is that the matrix c~(m) is rank deficient if the edge features are moving in a configuration such that
~ /3,~c11(rn) ~Ii
=
0
for all, m, with non-zero coefficients f3.~. For example, the algorithm fails if the image points used all happen to lie on a line in the image plane as such a line only defines a line on the visible surface from which it is impossible to infer the orientation of the surface [4]. A less trivial failure occurs if the image points lie on a smooth contour f(x, y) 0 that is a conic section [35]. Similar algorithms derived recently by Waxman and Wohn also fail in this case [37]. =
5. Conclusion Algorithms for interpreting the motion of edge features in an image sequence in terms of the position, orientation and motion of a planar visible surface facet have been developed and implemented on an SIMD processor array. Loosely speaking, these algorithms represent a first attempt at deriving a scheme for using visual motion to obtain 3D scene information at the level of Marr’s 2~Dsketch [38]. The scheme can be implemented very simply and efficiently on an SIMD processor array such as the DAP and it has been shown that, for the ideal 1 : 1 mapping of the image pixels onto the processor array, the algorithms run at speeds that are nearly compatible with real time video frame rates. An important problem remaining however, is segmenting the scene in order to be able to describe a situation in which there may be several objects moving against the background, or in which there may be objects of complex shape with several faces in view etc. A generalization of our procedure [32] which would achieve such a segmentation is currently under investigation.
280
B.F Buxton et al.
/
Structure-from - motion algorithms on an SIMD architecture
Acknowledgements
[14] H.C. Longuet-Higgins, Nature 293 (1981) 133. [15] H.C. Longuet-Higgins, Proc. Roy. Soc. London B223
It is a pleasure to thank Christopher LonguetHiggins and Guy Scott for many useful discussions and much advice. Particular thanks go to Professor Longuet-Higgins for timely preprints of his analysis of the interpretation of the visual motion due to a moving plane that stimulated our work on the 3D solutions to the aperture problem. This work has been carried out with the support of Procurement Executive, Ministry of Defence.
(1984) 165. [16] R.Y. Tsai and T.S. Huang, IEEE Trans. ASSP-29 (1981) 1147 R.Y.Tsai, T.S. Huang and W-L. Zhu, IEEE Trans. ASSP30 (1983) 525. R.Y. Tsai and T.S. Huang, IEEE Trans. ASSP-32 (1984) 213. [17] S.J. Maybank, in: Proc. ECAI 84, The Sixth Eur. Conf. on Artifical Intelligence (ECCAI, 1984) p. 641. [18] R. Paquin and E. Dubois, CVGIP 21(1983) 205.
References [1] J.M. Brady, Comput. Surv. 14 (1982) 3. [2] A.P. Reeves, CVGIP 25 (1984) 66. [31J.M. Brady, Artificial Intelligence 21 (1983) 271. [4] B.F. Buxton, H. Buxton, D.W. Murray and N.S. Williams, in: Proc. ECAI 84, The Sixth Eur. Conf. on Artificial Intelligence (ECCAI, 1984) p. 631. [5] S.F. Reddaway, Proc. 1st. Ann. Symp. on Computer Architecture (1973) p. 61—65. [6] D. Parkinson, Comput. Phys. Cominun. 28 (1983) 325. [7] D.J. Hunt and S.F. Reddaway, in: The Fifth Generation Computer Project, ed. 0G. Scarott (Pergamon Infotech, Maidenhead, 1982). [8] D.W. Lee, Phil. Trans. Roy. Soc. B290 (1980) 169. [9] H.C. Longuet-Higgins and K.F. Prazdny, Proc. Roy. Soc. B208 (1980) 385. [10] I. Hadam, G. Ishal and M. Gur, J. Opt. Soc. Am. 70 (1980) 60. [11] J.J. Koenderink and A.J. van Doom, Opt. Acta 22 (1975) 773. [12] AM. Waxman, in: Proc. Workshop on Computer Vision, Annapolis, 30 Apmil—2 May (IEEE Computer Soc. IEEE 84CH2014-9, 1984) p. 49. [13] S. Ullman, The Interpretation of Visual Motion (MIT Press, Cambridge, MA, 1979).
[19] [20] [21] [22]
H-H. Nagel, CVGIP 21(1983) 85. D. Marr and T. Poggio, Proc. Roy. Soc. B204 (1979) 301. W.E.L. Grimson, Phil. Trans. Roy. Soc. B292 (1981) 217. J.E.W. Mayhew and H.C. Longuet-Higgins, Nature 297 (1982) 376. [231 J.E.W. Mayhew, Perception 11 (1982) 387. [24] H.H. Baker, in: Proc. IEEE Intern. Conf. on Cybernetics in Society (1982). [25] B.K.P. Horn and B.G. Schunck, Artificial Intelligence 17 (1981) 185. [26] K.F. Prazdny, in: Proc. 7th Intern. Joint Conf. on Artificial Intelligence (1981) p. 698. [27] K.F. Prazdny, in: Proc. Conf. on Pattern Recognition and Image Processing, IEEE CH1761-6 (1982) p. 256. [28] B.F. Buxton and H. Buxton, Image and Vision Computing 2 (1984) 59. [29] B.F. Buxton and H. Buxton, GEC J. Res. 1 (1983) 184. [30] B.F. Buxton, H. Buxton and BK. Stephenson, lEE Proc. F 131 (1984) 593. [31] J.S. Wiejak, H. Buxton and B.F. Buxton, submitted to CVGIP. [32] G.L. Scott, private communication. [33] A.R. Bruss and B.K.P. Horn, CVGIP 21(1983) 3. [34] G.L. Scott, private communication. [35] B.F. Buxton, H. Buxton, D.W. Murray and N.S. Williams, submitted to GEC J. Res. [36] S. Ullman, MIT AIM-706 (1983). [37] A.M. Waxman and K. Wohn, CAR-TR-58, CS-TR-1394 (University of Maryland, 1984). [38] D. Marr, Vision (Freeman, San Francisco, 1982).
273
STRUCTURE-FROM-MOTION ALGORITHMS FOR COMPUTER VISION ON AN SIMD ARCHITECTURE B.F. BUXTON and D.W. MURRAY GEC Research Laboratories, Hirst Research Centre, East Lane, Wembley, UK
H. BUXTON and N.S. WILLIAMS Department of Computer Science and Statistics, Queen Mary College, Mile End Road, London El 4NS, UK
Algorithms for interpreting the motion of edge features in an image sequence in terms of the position, orientation and motion of a planar visible surface facet have been developed and implemented on an SIMD processor array. The underlying theory of the interpretation based on a simultaneous solution of the structure-from-motion and aperture problems is briefly reviewed and it is explained how the solution can be implemented very simply and efficiently on an SIMD processor array. Test calculations were carried out on the ICL DAP at Queen Mary College to illustrate the implementation and to show that, for an ideal 1 : 1 mapping of the image pixels onto the processor array, the algorithms run at speeds that are nearly compatible with real time video frame rates.
1. Introduction: computer vision Most of us are familiar with the idea that depth perception can be attained by stereo vision. Each eye (or sensor) sees a different view of the world with the result that the images formed are slightly out of alignment. From this shift or ‘disparity’ the depth can be computed with the help of a little geometry. However, it is not so widely appreciated that depth information can also be obtained if the sensor and the objects in view are moving relative to each other. One way to see this is to regard image snapshots taken at successive positions as stereo pairs. The two problems, both involving reconstruction of the depth dimension ‘missing’ in an ordinary image are obviously related. Moreover, in both cases it is the interpretation of the images in .terms of the three dimensional scene in view that generally distinguishes these problems as computer vision rather than belonging to the realm of image processing and pattern recognition. In image processing or pattern recognition the image is analysed as a two dimensional pattern per se, even when it is the projection of three dimensional objects in a three dimensional world [1].
The application of vector and parallel processors to image processing is now fairly commonplace, especially if SIMD or pipeline machines are used, and a number of more complicated architectures have been suggested for pattern recognition tasks (see ref. [21 for a recent review). Application of parallel computers in computer vision is not so widespread, but SIMD machines have been used for some low-level tasks and MIMD structures have been proposed for higher level vision work. Brady [3]has recently given a survey of parallelism in computer vision, including reference to some dedicated hardware implementations. In this paper, we want to describe the algorithms we have recently derived for obtaining 3D information from visual motion [4] and discuss their implementation on one of the simplest parallel architectures, an SIMD machine. The results of test calculations on the ICL DAP at Queen Mary College [5—7]will be used to illustrate the details of the implementation. The paper is organized as follows. We begin in the next section with a brief description of how optic flow (the technical term used to describe the motion of points in the image plane of the sensor)
OO1O-4655/85/$03.30 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
B.F Buxton et a!. / Structure-from - motion algorithms on an SIMD architecture
274
contains information about the depth and orientation of the visible surfaces and their motion relative to the sensor. We also explain why it is convenient to compute the motion of edge features in an image sequence rather than the whole of the optic flow field. As this only gives partial information about the flow field, we have to solve an aperture problem arising from the lack of information about movement of the edge features along their length in addition to reconstructing the missing depth information. The solution obtained recently [4] is outlined in section 3 and implementation of the algorithms derived is described in section 4. Particular attention is paid to the very simple way in which these algorithms map onto an SIMD processor array and the very high processing speeds that can be achieved. Geometrical conditions that cause the algorithms to fail and the numerical stability of the calculations are also briefly discussed. We end with a short conclusion and discussion that includes some comments on the direction of future work.
motion over the whole of the image plane have been discussed in a number of papers [8—12],and it has been shown that, in principle, the optic flow field contains enough information to allow one to use it to reconstruct the structure of rigid bodies moving relative to the sensor [13,14]. Obtaining the structure of a rigid body this way is referred to as the structure-from-motion problem. In our work, we consider a slight variant of the structure-from-motion problem in which we attempt to reconstruct visible surfaces from the visual motion. To show that the optic flow field contains information about the position and orientation of the visible surfaces, it is convenient to assume that the object points, R, lie on a planar surface facet, R N 1 2 3 — —
moving (instantaneously) with rectilinear velocity V and angular velocity ~Q with respect to the sensor. The optic flow field, r(r), given by 1=
V(rN)+r(~ V)(rN)/l+Q/\r +r[~.(QAr)]/l,
(2.4)
2. Optic flow As mentioned above, optic flow is the motion of points in the image plane due to motion of the corresponding object points relative to the sensor. The velocity of an image point r depends on the position, R (X, Y, Z), and velocity, R, of the object point according to (2.1) t 1R/Z+ IRZ/Z2, =
if we use the perspective transformation =
—
IR/Z
(2.2)
to model the geometry of the sensor (fig. 1). Details of the optic flow field, r(r), describing the
I
~J Fig. 1. The sensor geometry.
________>~
then depends on the surface and motion parameters N, V and £2. Since the terms involving N and V in (2.4) always appear as products of their components, we note that there is an inherent ambiguity in the optic flow field as it is the same for a small object close to, moving slowly, as for a large object, far away, moving quickly. Given the optic flowber at a small number of image points, (2.4) can solved for the surface and motion parameters [15—17]to within this inherent ambiguity (and a few other complications that we shall mention later in section 3). Thus one c,,an determine the direction of rectilinear motion V. the surface orientation, N, the angular velocity £2 and the the value of the products VN or v-N [15]. However, before these algorithms can be used, observed dE~)O th~ doing thisimage is to follow the motion of prominent image features from frame-to-frame, but this involves identifying prominent features in one frame and finding the corresponding features that they
B.F. Buxton et at
/
Structure-from - motion algorithms on an SIMD architecture
match in later frames. Algorithms of this kind involving a solution to the correspondence problem [18,19] are closely related to binocular stereo algorithm [20—24].There is however a simpler alternative approach, based on the motion constraint equa(ion [25—26], 3E ‘at —r E 2 ~ —
/
—
~‘
‘
/
‘
that avoids the correspondence problem. Since stands for the gradient operator (a/ax, a/ay) in the image plane, eq. (2.5) implies that the total derivative d E/d t is zero. Using the motion constraint equation is therefore equivalent to following the motion of a given image feature point which is usually assumed to be the same as obtaining the optic flow defined geometrically above, A fundamental difficulty with (2.5) however is that it can only yield the component of the optic flow parallel to ~E denoted here by r~as it is perpendicular to the edge contour. This lack of information about motion perpendicular to ç7E is known as the aperture problem, but the way it manifests itself depends on how (2.5) is used. Since the image irradiance values E ( r, t) are always corrupted by sensor noise, smoothing or filtering is required before (2.5) can be used. In addition, (2.5) will only give reliable estimates of i~where both aE/at and lyE I are large as is the case at moving edge features in the image [26—28].As described elsewhere [28,29], we therefore smooth the image irradiance values E(r, I) over both space and time with a Gaussian filter G(r, t) ~ 2~2±~2t2),~’2a2 and apply (2.5) at the edge feae_ tures defined by the zero-crossings of the visual signal [28,29], ~‘
S(r,
t)*E(r,
t)= _(~2+i~)G(r,
275
out the spatio-temporal convolution in (2.6), detecting the zero-crossings in S and estimating the vernier velocities v from (2.7), have already been described [28—30].We need only note here that as similar operations are to be carried out at every pixel in the image, these calculations map very simply onto an SIMD processor array such as the DAP. In addition, the masking facilities available on the DAP enable us to write simple uniform algorithms that only calculate the vernier velocity v at the edge features where it is required [30]. Linear least-squares interpolation enables us to calculate both v and the positions of the edge features to pub-pixel accuracy. As shown in fig. 2 for a toy calculation simulating a bright circular blob painted on a plane moving directly towards the sensor, the vernier velocities then accurately depict the edges of features on the visible surface. Other examples are given in the refs. [4,28—30]. The final point to note here is that these low-level calculations are extremely fast on the DAP, even when not fully optimized [28,30]. For the ideal 1: 1 mapping of the image pixels directly into the processor array, calculations of the kind illustrated in fig. 2 can be carried out in approximately 35 ms per frame, even though with the masks of radius 8 used (o 2%/L u I in the Gaussians), we are =
=
\\\
~ 1111 /f/ / / 1/,
.-
:
t).
(2.6) This yields an average value of the optic flow normal to the edge features, called the vernier velocity v [28]:
~
/
\“
~1, ~=~=
[(_as/at)(vs)] V
.
(2.7)
11,1 !•l~ ~ . Fig. 2. The vermer velocities at the edge of a brigbt circular
carrying
blobsensor. painted on the surface of a plane moving directly towards the
edge
Details of this calculation which involves
276
B. F. Buxton et a!.
/
Structure -from - motion algorithms on an SIMD architecture
effectively performing convolutions over 2000 pixels [30,31].
3. Interpreting edge motion The problem then is to interpret the edge motion by simultaneously solving the aperture and structure-from-motion problems. Eq. (2.4) relates the optic flow r to the motion and surface parameters V, £2 and N. To eliminate the unknown cornponent of the optic flow parallel to the image edge features, we take the dot product of (2.4) with the vernier v, leading to theV)(r.N)/l equation [4], 2 (v.velocity V)(r.N) +(vr)(z1 v +v.(S2Ar)+(v•r)[~(~2Ar)]/1, (3.1) =
for V, £2 and N. By writing the vectors in terms of their Cartesian components, V (V 1, V2, V3) etc., and using the summation convention, (3.1) may be rearranged in the form [4],
may be inverted for the parameters p,~leaving us to solve the eight non-linear simultaneous eqs. (3.4). The latter, however, is a problem that was encountered a little while ago by Longuet-Higgins in his work on the interpretation of the image of a moving planar surface [15], and we can take over the results of his analysis. Briefly, the solution proceeds by eliminating the unknowns that appear linearly in (3.4), the ~k, and constructing a cubic equation for the ‘odd-man-out’ amongst the remaining quantities, 173 N3. Once 173 N3 has been obtained from this cubic, we can solve for the remaining V ]VJ. and finally for the ~ Note, as mentioned 2, that because only 1’J~occur ininsection (3.4), the magnitudes of Vproducts and N cannot be fixed independently. In addition, the vectors V and N can be interchanged and their signs altered without affecting the value of the products 1’~N~giving four interpretations of the
=
[v1~~
=
—(v.r)~e3x1/x31p,1.
(3.2)
where we have introduced the parameters P,i (3.3) and defined x3 to be —l so that r (x1, x2, x3). As it stands, (3.2) looks like an equation for the nine parameters p~,but since the vernier velocity v always lies in the image plane, v3 is zero with the result that only the combinations ~ p33 and P22 p33 appear in (3.2). We therefore rewrite (3.2) as a set of linear simultaneous equations for the eight parameters, —
=
—
—
Pj1
=
~
—
~ijk~k
—
V3N36~1,
(3.4)
at each of the edge feature points, m, in matrix form [4]: 2(m). (3.5) ~ c~,(m)p11 v =
visual motion. Although only one of these would be consistent with the observed edge feature positions at later times, information immediately at hand can be used to reduce the number of possible interpretations. This is because we know that the surface facet described by N must be visible from each of the feature points R(m) used in the solution. R(m).N must therefore be positive leaving, at most, only two viable interpretations. If in addition, there are points m such that R(m) V> 0 and other points m’ for which R (rn’)’ V < 0, only one interpretation is consistent with all the data. Thus, Longuet-Higgins [15] states the theorem that there is a unique interpretation if there are visible points in both the forward and backward fields of view with respect to the direction of rectilinear motion V. As Longuet-Higgins [15] has shown that the middle root of the cubic equation gives the required value of V3N3 (the others, if distinct, give rise to complex flow fields) the solution looks straightforward, identity for cos 30 especially can be since used to thecalculate trigonometric V
3N3 directly. The only remaining problem is the solution of the simultaneous eqs. (3.5) for the parameters p,~,which still involves processing the data at
The matrix c~(m)is defined by: (3.6) With data at eight or more edge points, m, (3.6)
the edge points m. In principle, this too is simple enough: we either solve the set of eqs. (3.5), which may be written in matrix form for n edge features,
B.F Buxton et a!.
/ Structure -from
-
motion algorithms on an SJMD architecture
as
(ii) calculate the components of the coefficient 2
(3.7) directly in a least squares sense, or use the classical least squares solution, [c]5~g[p]8xj[v
I
277
Ti
]nxS,
r i.....r T21
v~cJgxsl.pJ
—
1C
V
J8X1~
.
/
Since the equations are not usually very well conditioned, the former would be preferable, but we have not found a convenient routine available to calculate the pseudo-inverse of c on the DAP (e.g. by Householder transformations as in the NAG routines) when the matrix is stored with its long dimension (n) mapped onto a DAP plane and its short dimension (8) mapped vertically. We have therefore used the classical least squares solution (3.8) as the data is arranged almost ideally for this calculation (section 4).
4. Implementation on the DAP Before describing the details of the implementalion of our algorithms on the DAP at Queen Mary College, it is convenient to summarize the computational steps involved in calculating the motion and surface parameters, V, £2, N and (NV). For simplicity we shall assume at this stage that all visible points lie on the same planar facet. In general, of course, this will not be true and we would like to use optic flow to segment the scene into objects moving against the background, visible surface facets in various orientations, etc. [28,29]. Current work [32] suggests that an iterative segmentation scheme based on the simple calculation we shall describe below can be devised. We comment further on this in the discussion (section 5), but note here that the present calculation, in which it is assumed that there is a single set of surface and motion parameters valid across the whole image, may be of direct relevance to passive navigation [33]. When only a single set of surface and motion parameters are required, we have to: (i) obtain the position of the edge features r( m) and calculate the vernier velocities v( m) describing the edge motion,
2], matrix c defined (iii) compute 1(m) [cTc]and [cTvby (3.6), (iv) solve the least squares eqs. (3.8) for the parameters p,~, (v) find the middle root, V 3 N3, of the cubic
equation, (vi) use V3N3 to compute the surface and motion parameters N, V, £2 and (VN), (vii) use the position of the edge features to resolve the interpretations of V and N as far as possible. Note that steps (i)—(iii) and (vii) involve calculations over the image space (x, y) whereas (iv), (v) and (vi) involve only the small space of the eight surface and motion parameters. In the former therefore, it is most convenient to map the image pixels into the processor array. This has already been described [30] for step (i) involving convolution of the images as in (2.6), detection of the zero-crossing edge features and computation of the edge positions r( m) and vernier velocities v( m). A sequence of small 64 X 64 images can be mapped directly onto the DAP, but for larger images we have to map groups of several neighbouring pixels onto the same processor in a ‘crinkled’ or ‘pyramidal’ mapping scheme [30]. In either case however, the vectors describing the details of the edge position and the vernier velocities end up stored vertically below the processors in the DAP mem-
64 x 64 processor array
16K of memory
Fig. 3. Sketch of the architecture of the 64 x 64 ICL DAP showing the connectivity of the processing elements each haying access to its own memory.
278
B.F Buxton et a!.
/
Structure-from - motion algorithms on an SIMD architecture
ory planes (fig. 3). For simplicity, therefore, in what follows we shall assume the ideal 1: 1 direct mapping of the image pixels onto the DAP. Only small changes are required for larger images. We should note however that the masking facilities on the DAP are used to isolate the edge pixels so that no explicit encoding of the edge positions is required [30]. The r(m) are used to record the position of the zero-crossings to subpixel accuracy [28—30]. With the positions r( m) and the vernier velocities v( m) stored vertically below each edge pixel, it is then very easy to calculate the components of the coefficient matrix c,~(m)from (3.6) and store these vertically below each processor as well (fig. 4), where they are ideally positioned for computing [cTc] and [cTv2] in step (iii). This is because I cTc i ~ (~m)\ c, / \ ~ 1 ~m, — —
m
and 2] [cTv
(4.1) cj~(m) v2 (rn),
=
ni
enabling us to carry out the multiplications in parallel using planar multiplication on the DAP, and to perform the sums, ~m’ over the edge pixels using the planar sum function. The result is a very efficient routine that can be programmed in four or five lines of DAP Fortran [6,7], made even faster by the fact that cTc is, of course, a symmetric matrix. A statistical weighting w(m) of each point to improve the least squares fit as suggested by Scott [4,34,35] can also be incorporated with only minor changes.
_.-‘I~~’~ Columns of c
Processor array rows,m,of c”
,_‘7]
The next three steps, (iv)—(vi) in which the surface and motion parameters are computed could be carried out on an attached processor (the host machine in the case of the DAP), although in our test computations the least squares eqs. (3.8), step (iv), were also solved on the DAP in order to exploit the parallelism in the Gauss—Jordan elimination [6]. Steps (v) and (vi) were carried out on the host, however, and the results transferred back to the DAP for step (vii). Approximate times per frame for the steps carried out on the DAP are summarized in table 1, assuming as mentioned above the ideal 1: 1 mapping of the image pixels onto the processor array. All calculations were carried out in REAL *4 arithmetic except for the spatio-temporal convolution involved in step (i) which is done in 32-bit integer arithmetic (and takes roughly 20 ms, only half the interframe time at normal video rates). Note that a little time is ‘wasted’ on converting data from the DAP storage mode to the host 2980 storagebemode and back out again, operations which would better carried in separate hardware in a dedicated vision system [30]. Timings for steps (v) and (vi) on the host ICL 2980 are not given because in a vision system these calculations could be carried out on an attached serial processor (say) whilst the heavily pressed SIMD processor array continues with other calculations on steps (i) and (iii) (say). This kind of
Table 1 Approximate times per frame for the steps carried out on the DAP Step Function Time Machine (ms) (i) (ii)
(iv) (iii)
~ k L
(v) (vi) Fig. 4. Storage of the matrix c, 1(m) in DAP memory showing how row indices, m, map into the processor array and how the columns q and k! (say) are stored vertically in the memory planes.
(vii)
calculate r(m) and v(m) set up c. Tcl and [cTv2l 3(m) solve for the p,, calculate[c DAP-2980 storage mode conversion solve the cubic compute surface and motion parameters 2980-DAP storage mode conversion find correct interpretation
36
DAP
354 8 7
DAP DAP DAP
HOST HOST 3 4
DAP
B.F Buxton et a!. / Structure-from - motion algorithms on an SIMD architecture
concurrency would be possible, of course, since in applications video data would be continuously arriving for new processing to begin. Similarly, although we solved the least square equations (step (iii)), on the DAP, this could also be hived off to an attached processor as it only involves the small motion and surface parameter space. However, this question requires more careful consideration for a complicated scene when there would (after segmentation) be several sets of motion and surface parameters to be considered. Finally, we should consider the numerical stability of our algorithms. Solving for the surface and motion parameters from the p,~appears to be quite stable unless two roots of the cUbic for V~ N3 are coincident or nearly so. According to Longuet-Higgins’ analysis [15] this happens when V and N are parallel or antiparallel and the sensor is moving directly towards or away from the plane. Given the pt,, steps (v) and (vi) are thus usually quite stable but we remarked above and in section 3 that the least squares fit for the p~ is not well-conditioned. There are two reasons for this. The first is that all 3D scene interpretations derived from structure-from-motion algorithms seem to depend to some extent on perspective effects that will be small unless the points in view encompass a large visual angle [36]. This can be seen directly for our algorithm by referringback to eq. (3.1) and assuming that the image co-ordinates x and y are small compared to the camera length 1. Since r (x, y, —1), (3.1) reduces in lowest order to =
= =
IN3 ( v1V1 + v2V2) lv1Q2 + lv2~21 —1v1(V1N3 + £22) —1v2( V2N3 ~ —
—
—
(4.2)
showing that only the parameters P13 and P23 can be reliably determined if the field of view is small, This shows that all but two of the singular values of the matrix c.~(m)will be small as may be verified directly. It is also in agreement with our previous results [4] on the stability of these structure-from-motion algorithms when the motion parameters are known a priori and we have to solve (3.1) for the surface normal N. To see this we note that if V and £2 are known, we can use the values of
p13=V1N3+Q2
and
p23=V2N3—f2~
(4.3)
279
to estimate N3 accurately, as found in ref. [4]. The second reason why (iv) can be unstable is that the matrix c~(m) is rank deficient if the edge features are moving in a configuration such that
~ /3,~c11(rn) ~Ii
=
0
for all, m, with non-zero coefficients f3.~. For example, the algorithm fails if the image points used all happen to lie on a line in the image plane as such a line only defines a line on the visible surface from which it is impossible to infer the orientation of the surface [4]. A less trivial failure occurs if the image points lie on a smooth contour f(x, y) 0 that is a conic section [35]. Similar algorithms derived recently by Waxman and Wohn also fail in this case [37]. =
5. Conclusion Algorithms for interpreting the motion of edge features in an image sequence in terms of the position, orientation and motion of a planar visible surface facet have been developed and implemented on an SIMD processor array. Loosely speaking, these algorithms represent a first attempt at deriving a scheme for using visual motion to obtain 3D scene information at the level of Marr’s 2~Dsketch [38]. The scheme can be implemented very simply and efficiently on an SIMD processor array such as the DAP and it has been shown that, for the ideal 1 : 1 mapping of the image pixels onto the processor array, the algorithms run at speeds that are nearly compatible with real time video frame rates. An important problem remaining however, is segmenting the scene in order to be able to describe a situation in which there may be several objects moving against the background, or in which there may be objects of complex shape with several faces in view etc. A generalization of our procedure [32] which would achieve such a segmentation is currently under investigation.
280
B.F Buxton et al.
/
Structure-from - motion algorithms on an SIMD architecture
Acknowledgements
[14] H.C. Longuet-Higgins, Nature 293 (1981) 133. [15] H.C. Longuet-Higgins, Proc. Roy. Soc. London B223
It is a pleasure to thank Christopher LonguetHiggins and Guy Scott for many useful discussions and much advice. Particular thanks go to Professor Longuet-Higgins for timely preprints of his analysis of the interpretation of the visual motion due to a moving plane that stimulated our work on the 3D solutions to the aperture problem. This work has been carried out with the support of Procurement Executive, Ministry of Defence.
(1984) 165. [16] R.Y. Tsai and T.S. Huang, IEEE Trans. ASSP-29 (1981) 1147 R.Y.Tsai, T.S. Huang and W-L. Zhu, IEEE Trans. ASSP30 (1983) 525. R.Y. Tsai and T.S. Huang, IEEE Trans. ASSP-32 (1984) 213. [17] S.J. Maybank, in: Proc. ECAI 84, The Sixth Eur. Conf. on Artifical Intelligence (ECCAI, 1984) p. 641. [18] R. Paquin and E. Dubois, CVGIP 21(1983) 205.
References [1] J.M. Brady, Comput. Surv. 14 (1982) 3. [2] A.P. Reeves, CVGIP 25 (1984) 66. [31J.M. Brady, Artificial Intelligence 21 (1983) 271. [4] B.F. Buxton, H. Buxton, D.W. Murray and N.S. Williams, in: Proc. ECAI 84, The Sixth Eur. Conf. on Artificial Intelligence (ECCAI, 1984) p. 631. [5] S.F. Reddaway, Proc. 1st. Ann. Symp. on Computer Architecture (1973) p. 61—65. [6] D. Parkinson, Comput. Phys. Cominun. 28 (1983) 325. [7] D.J. Hunt and S.F. Reddaway, in: The Fifth Generation Computer Project, ed. 0G. Scarott (Pergamon Infotech, Maidenhead, 1982). [8] D.W. Lee, Phil. Trans. Roy. Soc. B290 (1980) 169. [9] H.C. Longuet-Higgins and K.F. Prazdny, Proc. Roy. Soc. B208 (1980) 385. [10] I. Hadam, G. Ishal and M. Gur, J. Opt. Soc. Am. 70 (1980) 60. [11] J.J. Koenderink and A.J. van Doom, Opt. Acta 22 (1975) 773. [12] AM. Waxman, in: Proc. Workshop on Computer Vision, Annapolis, 30 Apmil—2 May (IEEE Computer Soc. IEEE 84CH2014-9, 1984) p. 49. [13] S. Ullman, The Interpretation of Visual Motion (MIT Press, Cambridge, MA, 1979).
[19] [20] [21] [22]
H-H. Nagel, CVGIP 21(1983) 85. D. Marr and T. Poggio, Proc. Roy. Soc. B204 (1979) 301. W.E.L. Grimson, Phil. Trans. Roy. Soc. B292 (1981) 217. J.E.W. Mayhew and H.C. Longuet-Higgins, Nature 297 (1982) 376. [231 J.E.W. Mayhew, Perception 11 (1982) 387. [24] H.H. Baker, in: Proc. IEEE Intern. Conf. on Cybernetics in Society (1982). [25] B.K.P. Horn and B.G. Schunck, Artificial Intelligence 17 (1981) 185. [26] K.F. Prazdny, in: Proc. 7th Intern. Joint Conf. on Artificial Intelligence (1981) p. 698. [27] K.F. Prazdny, in: Proc. Conf. on Pattern Recognition and Image Processing, IEEE CH1761-6 (1982) p. 256. [28] B.F. Buxton and H. Buxton, Image and Vision Computing 2 (1984) 59. [29] B.F. Buxton and H. Buxton, GEC J. Res. 1 (1983) 184. [30] B.F. Buxton, H. Buxton and BK. Stephenson, lEE Proc. F 131 (1984) 593. [31] J.S. Wiejak, H. Buxton and B.F. Buxton, submitted to CVGIP. [32] G.L. Scott, private communication. [33] A.R. Bruss and B.K.P. Horn, CVGIP 21(1983) 3. [34] G.L. Scott, private communication. [35] B.F. Buxton, H. Buxton, D.W. Murray and N.S. Williams, submitted to GEC J. Res. [36] S. Ullman, MIT AIM-706 (1983). [37] A.M. Waxman and K. Wohn, CAR-TR-58, CS-TR-1394 (University of Maryland, 1984). [38] D. Marr, Vision (Freeman, San Francisco, 1982).