Time-to-contact estimation device and method for estimating time to contact

Abstract

A time-to-contact estimation device for estimating a time to contact of a moving body to an object in surroundings of the moving body, includes: an image sensor provided on the moving body, and continuously capturing time-varying images of the surroundings of the moving body; a brightness change calculation unit calculating brightness changes in pixels between the time-varying images captured one after another; and a time-to-contact calculation unit calculating optimized brightness gradients based on the brightness changes calculated by the brightness change calculation unit, and calculating the time to contact of the moving body to the object based on the optimized brightness gradients.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic diagram illustrating a method for estimating the time to contact using a brightness gradient.

FIG. 2 is a schematic diagram illustrating definition of coordinate systems.

FIG. 3 is a schematic diagram illustrating conception of optical flow.

FIG. 4 is a schematic diagram illustrating conception of brightness acquisition.

FIG. 5 is a schematic diagram illustrating a method for calculating derivatives of brightness.

FIG. 6 is a schematic diagram illustrating expression of the derivatives of brightness by a cubic pixel.

FIG. 7 is a schematic diagram illustrating the relationship between a planer surface and a camera in Case 1 in which the degrees of freedom of the planar surface and the camera are constrained.

FIG. 8 is a schematic diagram illustrating the relationship between a planer surface and a camera in Case 2 in which the degree of freedom of the planar surface is constrained.

FIG. 9 is a schematic diagram illustrating the relationship between a planer surface and a camera in Case 3 in which the degree of freedom of the camera is constrained.

FIG. 10 is a schematic diagram illustrating the relationship between a planer surface and a camera in Case 4 in which none of the degrees of freedom of the planar surface and the camera are constrained.

FIG. 11 is a schematic diagram illustrating an estimation flow for the time to contact in Case 1.

FIG. 12 is a schematic diagram illustrating an estimation flow for the time to contact in Case 2.

FIG. 13 is a schematic diagram illustrating an estimation flow for the time to contact in Case 3.

FIG. 14 is a schematic diagram illustrating an estimation flow for the time to contact according to Solution 1 in Case 4.

FIG. 15 is a schematic diagram illustrating a calculation flow for the time to contact according to Solution 2 in Case 4.

FIG. 16 is a schematic diagram illustrating estimation of the time to contact using a feature tracking method.

FIG. 17 is a schematic diagram illustrating problems in the feature tracking method.

DETAILED DESCRIPTION OF THE INVENTION

Hereinafter preferred embodiments of the invention will be described with reference to the appended drawings.

FIG. 1 is a schematic diagram illustrating a method for estimating the time to contact using a brightness gradient. FIG. 1 illustrates the case in which a subject vehicle 10 (a moving body) is approaching another vehicle 20 (an object) traveling across the road on which subject vehicle 10 proceeds.

FIG. 2 is a schematic diagram illustrating definition of coordinate systems, in which the position p of the object is expressed by coordinates (x,y) in an imaging system oriented coordinate system, and the position P of the object is expressed by coordinates (X,Y,Z) in a sensor oriented coordinate system.

FIG. 3 is a schematic diagram illustrating conception of optical flow. As shown in FIG. 3, the optical flow is defined by a motion vector drawn between consecutive images of an object, and is used to detect or measure motion of an object. Practically, the optical flow is calculated as a motion vector drawn between two points under assumption that a first point in a first image appears as a second point having the same brightness as the first point and shifted from the first point by a small amount in a subsequent second image. As shown in FIG. 4, the brightness at time t of a point P(x,y) in an image is expressed by E(y,y,t). Assuming that the point P moves by (Δx,Δy) from time t to time t+Δt while maintaining brightness thereof, following Equation (1) is satisfied.

E(x,y,t)=E(x+Δx,y+Δy,t+Δt)   (1)

The Taylor's deployment is applied to the right side of Equation (1) and the second or higher order terms are neglected to obtain following Equation (2).

$\hspace{1em}\begin{array}{cc}\begin{array}{c}E\left(x,y,t\right)\cong E\left(x,y,t\right)+\frac{\partial E\left(x,y,t\right)}{\partial x}\Delta x+\frac{\partial E\left(x,y,t\right)}{\partial y}\Delta y+\\ \frac{\partial E\left(x,y,t\right)}{\partial t}\Delta t\end{array}& \left(2\right)\end{array}$

By dividing both sides by At to obtain following Equation (3).

$\begin{array}{cc}\frac{\partial E\left(x,y,t\right)}{\partial x}\frac{\Delta x}{\Delta t}+\frac{\partial E\left(x,y,t\right)}{\partial y}\frac{\Delta y}{\Delta t}+\frac{\partial E\left(x,y,t\right)}{\partial t}=0& \left(3\right)\end{array}$

The motion vector in the image is defined by following Equation (4);

$\begin{array}{cc}\left(u,v\right)=\left(\frac{x}{t},\frac{y}{t}\right)& \left(4\right)\end{array}$

the spatial gradient of the brightness is defined by following Equation (5); and

$\begin{array}{cc}\left(E_x,E_y\right)=\left(\frac{E\left(x,y,t\right)}{\mathrm{dx}},\frac{E\left(x,y,t\right)}{\mathrm{dy}}\right)& \left(5\right)\end{array}$

the temporal gradient of the brightness is defined by following Equation (6).

$\begin{array}{cc}{E}_t=\frac{E\left(x,y,t\right)}{\mathrm{dt}}& \left(6\right)\end{array}$

Equation (3) is converted by using Equations (4), (5), and (6) to obtain following Equation (7).

E _x ·u+E _y ·v+E _t=0   (7)

Equation (7) is a brightness gradient constraint equation. Parameters u and v that satisfy Equation (7) define the optical flow, i.e., the motion vector to be obtained.

As shown in FIG. 5, practically, E_x, E_y, and E_t are calculated as follows by using the brightness of pixels adjacent to each other in the image at time t and at time t+Δt, where i and j are arbitrary coordinates in the image.

E _x=(E(i+1,j,t)−E(i,j,t))+(E(i+1,j,t+Δt)−E(i,j,t+Δt))+(E(i+1,j+1,t)−E(i,j+1,t))+(E(i+1,j+1,t+Δt)−E(i,j+1,t+Δt))   (8)

E _y=(E(i,j+1,t)−E(i,j,t))+(E(i,j+1,t+Δt)−E(i,j,t+Δt))+(E(i+1,j+1,t)−E(i+1,j,t))+(E(i+1,j+1,t+Δt)−E(i+1,j,t+Δt))   (9)

E _t=(E(i,j,t+Δt)−E(i,j,t))+(E(i,j+1,t+Δt)−E(i,j+1,t))+(E(i+1,j,t+Δt)−E(i+1,j,t))+(E(i+1,j+1,t+Δt)−E(i+1,j+1,t))   (10)

E_x, E_y, and E_t can also be expressed by using conception of a cubic pixel as shown in FIG. 6 (refer to Horn & Schunck 1981).

In the present invention, in order to accurately estimate the time to contact, the positional relationship between the object (a planar surface 1 in the figures) and the image sensor (a camera 2 in the figures) is classified into four cases shown in FIGS. 7 to 10.

Note that, in the following description, the term “the degree of freedom of the camera 2” is defined by whether or not the camera 2 moves along the optical axis 2a thereof, i.e., the degree of freedom of the camera 2 is deemed to be constrained when the camera 2 moves along the optical axis 2a, and the degree of freedom of the camera 2 is deemed not to be constrained when the camera 2 moves not necessarily along the optical axis 2a. Furthermore, the term “the degree of freedom of the planar surface 1” is defined by whether or not the planar surface 1 is oriented perpendicular to the optical axis 2a, i.e., the degree of freedom of the planar surface 1 is deemed to be constrained when the planar surface 1 is oriented perpendicular to the optical axis 2a, and the degree of freedom of the planar surface 1 is deemed not to be constrained when the planar surface 1 is not necessarily oriented perpendicular to the optical axis 2a.

Case 1

In Case 1, as shown in FIG. 7, a translational motion of the camera 2 along a translation path 3, which is in parallel with the optical axis 2a of the camera 2 in this case, toward the planar surface 1 fixed perpendicular to the optical axis 2a (i.e., the normal 1a of the planar surface 1 is in parallel with the optical axis 2a and the translation path 3) is discussed. In this case, because the degrees of freedom of the planar surface 1 and the camera 2 are constrained, i.e., the camera 2 moves along the optical axis 2a thereof without rotational motion, the brightness gradient constraint equation is the simplest form. In addition, by assuming that the object to be considered is a planar surface that is oriented perpendicular to the optical axis 2a of the camera 2, motion vector (u,v) in the image (u=dx/dt, v=dy/dy) can be simply expressed as follows,

u=x(W/Z), v=y(W/Z)   C1-(1)

and the time to contact TTC is expressed as follows,

TTC=Z/W   C1-(2)

where W is speed of the camera 2 along the Z axis shown in FIG. 2, and Z is the distance to the planar surface 1.

By referring to above Equation (7), the brightness gradient constraint equation shown below is satisfied,

rE _r C+E _t=0   C1-(3)

where C(=W/Z) is the inverse of the TTC.

In addition, following Equations C1-(4) are satisfied, where r is the radial distance from the principal point.

rE _r =xE _x +yE _y , r=√{square root over (x²+y²)}  C1-(4)

The unknown parameter C can be found based on Equation C1-(3) using the least squares method. By squaring the both sides of Equation C1-(3) and taking the sum, following summed-up term C1-(5) is obtained.

Σ(rE_rC+E_t)²   C1-(5)

Summed-up term C1-(5) is over all pixels of a region of interest (ROI), or possibly over the whole image. To find the best fit value of C, summed-up term C1-(5) is differentiated with respect to C, and the result is set to be zero.

2Σ(rE_rC+E_t)rE_r=0   C1-(6)

Equations C1-(6) can be converted into following Equation C1-(7).

CΣ(rE_r)²=−ΣrE_rE_t   C₁-(7)

The unknown C, which is the inverse of the TTC, can be obtained by the following Equation C1-(8)

C=−ΣrE _r E _t/Σ(rE_r)²   C1-(8)

With reference to FIG. 11, the calculation steps of determining the TTC in Case 1 are as follows.

• (1) Using two (or more) frames from an image sequence, calculate the derivatives E_x, E_y, and E_t at every pixel in the region of interest from the brightness at the pixel in one image and the brightness at neighboring pixels in that image and the one before it (step C1-S1).
• (2) At each pixel, compute the term G=rE_r=xE_x+yE_y, from that compute the terms GE_t and G² (step C1-S2), and add these two terms, respectively.
• (3) Estimate the unknown C using Equation C1-(8) (step C1-S3).
• (4) Calculate the TTC as the inverse of C (step C1-S4).

The computation of the spatial and temporal derivatives preferentially uses a 2×2×2 cube of pixel brightness values (Horn & Schunck 1981). However, other methods of estimating derivatives can be used also, such as central differences, including methods that depend on more than two frames of the image sequence (although this may increase latency).

No significant amount of storage for intermediate results is needed, since the derivatives computed in step (1) are immediately used in step (2), and the terms computed in step (2) can be added into the accumulators in the immediately following step and so need not be retained.

Contributions to the accumulated totals can be weighted, so that, for example, more trustworthy information contributes more heavily to the final result. If this is done, the same weight factor should be applied to contributions to both accumulators. Information from an area of the image that has no significant texture, for example, may not provide a useful contribution to the accumulated totals. Also, areas in which specular highlights are observed may be suppressed in order to reduce error introduced by specular reflections.

It should be clear that the computation is very straightforward, with a number of arithmetic steps proportional to the number of pixels in the region of interest, that no search or iteration is involved and there is no need to solve complex equations.

Further, no detection of image “features” such as edges or grey-level corners is required. This makes it possible to implement the method on a small microprocessor (such as an Analog Devices BlackFin microprocessor) or even in special purpose circuitry such as FPGA, PLD or ASIC. In fact, the circuitry for performing the computation could even be integrated with the imaging sensor on the same. The same observations apply to the other cases studied below.

Case 2

In Case 2, as shown in FIG. 8, a translational motion of the camera 2 along a translation path 3, which is arbitrarily determined in this case, toward the planar surface 1 fixed perpendicular to the optical axis 2a (i.e., the normal 1a of the planar surface 1 is in parallel with the optical axis 2a) is discussed. In this case, the degree of freedom of the planar surface 1 is constrained, and the degree of freedom of the camera 2 is not constrained, i.e., the camera 2 moves not necessarily along the optical axis 2a, or not necessarily along a line in parallel with the normal 1a of the planar surface 1. U, V, and W shown in FIG. 2 respectively indicate X, Y, and Z component of velocity in the sensor oriented coordinate system. With regard to these terms, following Equation C2-(1) (i.e., brightness gradient constraint equation) is satisfied, where f is focal length of the camera 2.

u=−f·U/Z+x·W/Z, v=−f·V/Z+y·W/Z   C2-(1)

It should be clear that Case 1 discussed above is simply a special case of Case 2 with U=V=0.

Here, the following brightness gradient constraint equation is satisfied, where A=f·U/Z, B=f·V/Z, and C=W/Z (C is the inverse of the TTC as in above Case 1).

−AE _x −BE _y +rE _r C+E _t=0   C2-(2)

In addition, as in above Case 1, following Equation C2-(3) is satisfied.

G=rE _r =xE _x +yE _y , r=√{square root over (x²+y²)}  C2-(3)

The unknown parameters A, B, and C can be found based using the least squares method, where summed-up term is as follows.

Σ(−AE_x−BE_y+rE_rC+E_t)²   C2-(4)

As in Case 1, summed-up term C2-(4) is over all pixels of a region of interest (ROI), or possibly over the whole image. To find the best fit values of A, B, and C, summed-up term C2-(4) is differentiated with respect to A, B, and C, and the three results are set to be zero.

Σ(−AE_x−BE_y+rE_rC+E_t)E_x=0   C2-(5)

Σ(−AE_x−BE_y+rE_rC+E_t)E_y=0   C2-(6)

Σ(−AE_x−BE_y+rE_rC+E_t)rE_r=0   C2-(7)

These three linear equations containing the three unknown parameters A, B, and C can be solved easily using known methods such as Gaussian elimination.

With reference to FIG. 12, the calculation steps of determining the TTC in Case 2 are as follows.

• (1) Using two (or more) frames from an image sequence, calculate the derivatives E_x, E_y, and E_t at every pixel in the region of interest from the brightness at the pixel in one image and the brightness at neighboring pixels in that image and the one before it (step C2-S1).
• (2) At each pixel, compute the term G=rE_r=xE_x+yE_y(step C2-S2).
• (3) Estimate the unknown parameters A, B, and C based on Equations C2-(5), C2-(6), and C2-(7) (step C2-S3).
• (4) Calculate the TTC as the inverse of C (step C2-S4).

Note that the six distinct quantities in the symmetric 3×3 coefficient matrix, and the three quantities in the right hand side vector are all sums of products of various combinations of brightness derivatives Ex, Ey, and Et and image coordinates x and y.

The time to contact is the inverse of C=(W/Z). The focal length f (principal distance) need not be known to compute the time to contact. As noted, Case 1 discussed above is a special case of case 2, with A=B=0. In that special case, only the third of the three equation above needs to be considered, and two of its terms drop out.

In addition to the time to contact, if desired, one can calculate the direction of translational motion by noting that

U/W=(A/C)/f and V/W=(B/C)/f.

The magnitude of the translational motion, however, cannot be determined in this fashion due to the scale factor ambiguity.

The results can be computed using essentially the same sequence of eight steps as in Case 1 above, except that now 9 (6 plus 3) accumulators are needed instead of two, and nine terms per picture cell need to be computed instead of just two.

If desired, the stability of the solution can be analyzed by determining the eigenvalues and eigenvectors of the symmetric 3×3 coefficient matrix above. The result is well constrained if none of the eigenvalues are particularly small. Conversely, the accuracy can be expected to be low if at least one of the eigenvalues is small.

Case 3

In Case 3, as shown in FIG. 9, a translational motion of the camera 2 along a translation path 3, which is in parallel with the optical axis 2a of the camera 2 in this case, toward the planar surface 1 of an arbitrary orientation (i.e., the normal 1a of the planar surface 1 is not necessarily in parallel with the optical axis 2a) is discussed. In this case, the degree of freedom of the camera 2 is constrained, and the degree of freedom of the planar surface 1 is not constrained.

When p and q are defined to be the slopes in the X and Y direction of the planar surface 1 measured in the imaging system coordinate system, the planar surface 1 is expressed by the following Equation C3-(1).

Z=Z0+pX+qY   C3-(1)

X=Z·x/f and Y=Z·y/f are substituting into Equation C3-(1) to obtain following Equation C3-(2).

Z(1−p(x/f)−q(y/f))=Z0   C3-(2)

By using the relationship among brightness, these motion vector constraint equations can be expressed by following Equation C3-(3).

rE _r(1−p·x/f−q·y/f)(W/Z0)+E_t=0   C3-(3)

Equation C3-(3) is converted to obtain following Equation C3-(4), where C=W/Z0, P=(p/f)·C, and Q=(q/f)·C.

rE _r(C−Px−Qy)+E_t=0   C3-(4)

It should be clear that Case 1 discussed above is simply a special case of Case 3 with P=Q=0.

As in Cases 1 and 2, the unknown parameters P, Q, and C can be found based using the least squares method, where summed-up term is as follows.

Σ(rE_r(C−Px−Qy)+E_t)²   C3-(5)

As in Cases 1 and 2, summed-up term C3-(5) is over all pixels of a region of interest (ROI), or possibly over the whole image. To find the best fit values of P, Q, and C, summed-up term C3-(5) is differentiated with respect to A, B, and C, and the three results are set to be zero.

Σ(rE_r(C−Px−Qy)+E_t)rE_rx=0   C3-(6)

Σ(rE_r(C−Px−Qy)+E_t)rE_ry=0   C3-(7)

Σ(rE_r(C−Px−Qy)+E_t)rE_r=0   C3-(8)

These three linear equations containing the three unknown parameters P, Q, and C can be solved easily using known methods such as Gaussian elimination.

With reference to FIG. 13, the calculation steps of determining the TTC in Case 3 are as follows.

• (1) Using two (or more) frames from an image sequence, calculate the derivatives E_x, E_y, and E_t at every pixel in the region of interest from the brightness at the pixel in one image and the brightness at neighboring pixels in that image and the one before it (step C3-S1).
• (2) At each pixel, compute the term G=rE_r=xE_x+yE_y (step C3-S2).
• (3) Estimate the unknown parameters P, Q, and C based on Equations C3-(6), C3-(7), and C3-(8) (step C3-S3).
• (4) Calculate the TTC as the inverse of C (step C3-S4).

Note that the nine quantities in the symmetric 3×3 coefficient matrix, and the three quantities in the right hand side vector are all sums of products of brightness derivatives E_x, E_y, and E_t and image coordinates x and y.

The time to contact is the inverse of C=(W/Z0). The focal length f (principal distance) need not be known to compute the time to contact. It should be noted again that Case 1 discussed above is a special case of this one with P=Q=0. In this case, only the third of the three equations above need be considered, and two of its terms drop out.

If desired, one can also calculate the orientation of the planar surface by noting that

p=f(P/C) and q=f(Q/C).

The results can be computed using essentially the same sequence of eight steps as in case 1 above, except that now 9 (6 plus 3) accumulators are needed instead of two, and nine terms per picture cell need to be computed instead of just two.

If desired, the stability of the solution can be analyzed by determining the eigenvalues and eigenvectors of the symmetric 3×3 coefficient matrix above. The result is well constrained if none of the eigenvalues are particularly small. Conversely, the accuracy can be expected to be low if at least one of the eigenvalues is small.

The general approach should be clear from the above three cases worked out in detail above. Below described is an additional extension to more general case.

Case 4

In Case 4, as shown in FIG. 10, a translational motion of the camera 2 along a translation path 3, which is arbitrarily determined in this case, toward the planar surface 1 of an arbitrary orientation (i.e., the normal 1a of the planar surface 1 is not necessarily in parallel with the optical axis 2a) is discussed. In this case, none of the degrees of freedom of the camera 2 and the planar surface 1 are constrained.

Combining the equations for the motion field from Case 2 and the equation for the planar surface from Case 3 in the brightness change constraint equation, following Equation C4-(1) is obtained.

(A′E_x+B′E_y−rE_r)(−Px−Qy+C)+E_t=0   C4-(1)

where A′=f·(U/W),B′=f·(V/W),P=(p/f)·(W/Z0),Q=(q/f)·(W/Z0 )

As before, C=(W/Z0) is the inverse of the time to contact. This can also be written in the alternate form below.

(AE_x+BE_y−C(rE_r))(−P′x−Q′y+1)+E_t=0   C4-(2)

where A=f·(U/Z0),B=f·(V/Z0), P′=p/f,Q′=q/f

Note that Case 1 is just a special Case of case 4, with A=B=U=V=0, and correspondingly, Cases 2 and 3 are special cases of this general case with A=B=0 and U=V=0, respectively.

The least-squares method is applied to find the five unknown parameters A′, B′, C, P and Q, where summed-up term is as follows.

Σ((A′E_x+BE_y−rE_r)(−Px−Qy+C)+E_t)²   C4-(3)

Similarly, the least-squares method is applied to find the five unknown parameters A, B, C, P′ and Q′, where summed-up term is as follows.

Σ((AE_x+BE_y−C(rE_r))(−P′x−Q′y+1)+E_t)²   C4-(4)

Summed-up terms C4-(3) and C4-(4) are over all pixels of a region of interest (ROI), or possibly over the whole image. Note that A, B, P, Q are related to A′, B′, P′ and Q′ by the following equations.

A=A′C, B=B′C, P=P′C, Q=Q′C

To find the best fit values of the five unknown parameters we can differentiate either of the two sums with respect to the five parameters and set the results equal to zero. This leads to five equations in five unknowns. All of the coefficients of these equations are sums of products of image brightness derivatives and image coordinates, as before. Unlike Cases 1 to 3 described above, however, the equations are no longer linear and so not quite so easy to solve, and two solutions will be discussed below.

Solution 1

In Solution 1, an iteration process is applied to find unknown parameters defining the time to contact (TTC). First, parameters P′ and Q′ are treated as constants, then the equations derived from summed-up term C4-(4) are linear in the remaining unknown parameters A, B, and C, and known linear solution methods can be applied to obtain the best fit values of A, B, and C. Next, parameters A′ and B′ are treated as constants, then the equations derived from summed-up term C4-(3) are linear in the remaining unknown parameters P, Q, and C, and known linear solution methods can be applied to obtain the best fit values of P, Q, and C. By executing an iteration process for these calculations, approximation of the TTC can be obtained.

The specific procedure is now explained. First, parameters P′ and Q′ are treated as constants. To find the best fit values of A, B, and C, summed-up term C4-(4) is differentiated with respect to A, B, and C, and the three results are set to be zero.

AΣE _x ² F ² +BΣE _x ² E _y F ² −CΣrE _r E _x F ² =−ΣE _x E _t F   C4-(5)

AΣE _x ² F ² +BΣE _x ² E _y F ² −CΣrE _r E _x F ² =−ΣE _y E _t F   C4-(6)

AΣE _x ² F ² +BΣE _x ² E _y F ² −CΣrE _r E _x F ² =−ΣrE _r E _t F   C4-(7)

where F=1−(p/f)x−(q/f)y.

Next, parameters A′ and B′ are treated as constants. To find the best fit values of A, B, and C, summed-up term C4-(3) is differentiated with respect to P, Q, and C, and the three results are set to be zero.

PΣD ² x ² +QΣD ² xy−CΣD ² x=ΣD xE _t   C4-(8)

PΣD ² xy+QΣD ² y ² −CΣD ² y=ΣD yE _t   C4-(9)

PΣD ² x+QΣD ² y−CΣD ² =ΣDE _t   C4-(10)

where D=A′E_x+B′E_y−rE_r.

Thus, given an initial guess, one can alternately solve for A, B, and C, assuming P and Q are fixed, and then, using the new estimates of A, B, solve for P, Q and C, assuming that A and B are fixed. A few iterations of this pair of steps typically yields a close enough approximation to the exact solution.

With reference to FIG. 14, the calculation steps of estimating the TTC in Solution 1 are as follows.

• (1) Using two (or more) frames from an image sequence, calculate the derivatives (brightness gradients) E_x, E_y, and E_t at every pixel in the region of interest from the brightness at the pixel in one image and the brightness at neighboring pixels in that image and the one before it (step C4-1-S1).
• (2) At each pixel, compute the term G=rE_r=xE_x+yE_y (step C4-1-S2).
• (3) Calculate initial values (step C4-1-S3).
• (4) Estimate the unknown parameters P, Q, and C (step C4-1-S4).
• (5) Determine whether the solutions have been converged (step C4-1-S5). If the solutions have been sufficiently converged, the processing is terminated, and if not, the following steps are carried out.
• (6) Calculate p and q (step C4-1-S6).
• (7) Calculate F(=1−p/f·x−p/f·y) (step C4-1-S7).
• (8) Estimate the unknown parameters A, B, and C (step C4-1-S8), then return to step (3).

Solution 2

By using Solution 2, an approximation to the TTC can be obtained.

First, following terms C4-(11) and C4-(11) are derived from Equations C4-(1) and C4-(2).

(E_x E_y −rE_r)(A B C)^T(−P′ −Q′ 1)(x y 1)^T+E_t   C4-(11)

(E_x E_y −rE_r)(A′ B′ C)^T(−P −Q C)(x y 1)^T+E_t   C4-(12)

Suffix T denotes the transpose, which turns a row vector into a column vector (and vice versa). These expressions in turn can be rewritten using the product of a vector, a 3×3 matrix M, and another vector:

(E_x E_y −rE_r)M(x y 1)^T+E_t   C4-(13)

where M is the dyadic product of the vectors (A B C)^T and (−P′ −Q′ 1)^T, or equivalently of (A′ B′ 1)^T and (−P −Q C)^T, namely:

$\begin{array}{cc}M=\left(\begin{array}{ccc}-\mathrm{AP}& -\mathrm{AQ}& A\\ -\mathrm{BP}& -\mathrm{BQ}& B\\ -\mathrm{CP}& -\mathrm{CQ}& C\end{array}\right)=\left(\begin{array}{ccc}-A^{\prime }P& -A^{\prime }Q& A^{\prime }C\\ -B^{\prime }P& -B^{\prime }Q& B^{\prime }C\\ -P& -Q& C\end{array}\right)& \mathrm{C4}\text{-}\left(14\right)\end{array}$

A linear least squares method may now be used to find the nine unknown coefficients of the 3×3 matrix M, if these parameters are treated as independent variables for the moment. The term above can be written in the form (e^Tm+E_t), and so the sum to be minimized is

Σ(e^Tm+E_t)²   C4-(15)

where m is a vector obtained by “flattening” the matrix M, that is, by lining up the nine coefficients of the matrix, and e is a vector obtained by flattening the dyadic product of the vectors (−E_x −E_y (rE_r))^T and (x y 1)^T.

m=Flatten(((A B C)^T(−P′ −Q′ 1))   C4-(16)

m=Flatten(((A′ B′ 1)^T(−P −Q C))   C4-(17)

e=Flatten(((E_x E_y −rE_r)^T(x y 1))   C4-(18)

From Equations C4-(16), C4-(17), and C4-(18), following Equations C4-(19), C4-(20), and C4-(21) are obtained.

m=(−AP′ −AQ′ A −BP′ −BQ′ B −CP′ −CQ′ C)   C4-(19)

m=(−A′P −A′Q A′C −B′P −B′Q B′C −P −Q C)   C4-(20)

e=(xE_x yE_x E_x xE_y yE_y E_y −xrE_r −yrE_r −rE_r)   C4-(21)

By differentiating term C4-(15) with respect to the vector m, and setting the result equal to zero, following Equations C4-(22) and C4-(23) are obtained.

2Σ(e^Tm+E_t)e=0   C4-(22)

2Σ(ee^T)m=−ΣE_te   C4-(23)

This is a set of 9 linear equations in the nine coefficients of the matrix which can be solved using known methods, such as Gaussian elimination. The coefficients of the symmetric 9×9 matrix and the vector on the right hand side of the equation are all sums of products of derivatives of brightness and image coordinates. These terms can be accumulated as before, except that now there are 54 (45 plus 9) accumulators, and 54 terms per picture cell need to be computed for addition into these accumulators.

Finally, the five unknown parameters A, B, C, P′, and Q′ or equivalently A′, B′, C, P and Q can be recovered using singular value decomposition of the matrix M.

Alternatively, in this special case, they can be found by noting that each column of the matrix should be proportional to (A, B, C), and each row of the matrix should be proportional to (−P, −Q, C). In fact, the four unknown parameters A, B, P, and Q can be read off from the last row and the last column of the matrix M, reconstituted from the flattened version, the vector m.

However, the 3×3 matrix resulting from this computation may have an arbitrary factor lambda (λ) added to each of the three diagonal terms since

−xE _x −yE _y+(rE_r)=0.

The correct diagonal can be determined using singular value decomposition, or in this particular case more easily from relationships with the off diagonal terms:

M11=M12·M31/M32=M21·M13/M23   C4-(24)

M22=M12·M23/M13=M21·M32/M31   C4-(25)

M33=M32·M13/M12=M23·M31/M21   C4-(26)

Equation C4-(26) yields C, which is the inverse of the TTC. For purposes of determining the time to contact, only that bottom right coefficient of the matrix (Equation C4-(26)) actually is needed.

Solution 2 does not solve exactly the same least squares problem as that solved by Solution 1 described above, but does produce a result in convenient closed form. In the absence of measurement noise, the results will be the same. With noise in the measurement of image brightness, and hence in E_x, E_y and E_t, the result of Solution 2 will be an approximation to the original solution. This is because the nine elements of the matrix M are not independent, but depend on the five unknown parameters. This approximate solution may be sufficient enough to be used directly or can be used as an initial guess for the non-linear numerical solution of the original least-squares problem using iterative methods (i.e., for Solution 1).

As before, the time to contact is the inverse of the parameter C. This calculation can be performed without knowing the principal distance f.

If desired, the direction of translational motion, given by

U/W=(A/C)/f   C4-(27)

V/W=(B/C)/f   C4-(28)

can also be calculated, as can the orientation of the surface specified by

p=f(P/C) and q=f(Q/C).

With reference to FIG. 15, the calculation steps of estimating the TTC in Solution 2 are as follows.

• (1) Using two (or more) frames from an image sequence, calculate the derivatives (brightness gradients) E_x, E_y, and E_t at every pixel in the region of interest from the brightness at the pixel in one image and the brightness at neighboring pixels in that image and the one before it (step C4-2-S1).
• (2) Define vector e (step C4-2-S2).
• (3) The least-squares method is applied (step C4-2-S3).
• (4) Calculate diagonal elements of matrix M (step C4-2-S4).
• (5) Calculate each of the elements of vector m (step C4-2-S5).
• (6) Estimate the TTC (step C4-2-S6).

The specific embodiments described above are exemplary only and it will be apparent to the person of ordinary skill in the art that numerous alternatives exist as to the exact details of how the time to contact is computed using derivatives of brightness. Furthermore, other combinations of translational and rotational motion and other surface shapes can be treated using this method.

If there are multiple moving objects, the images can be segmented into regions of interest corresponding to images of individual objects.

The time to contact is often considered to be the time until the center of projection of the sensor intersects the surface of the object being viewed. If the sensor is mounted on an extended object of known shape, then the time when that shape first intersects with the object being viewed can be determined.

In the computation of the time to contact it is typically assumed that the motion parameters do not change. If the sensor or the object on which it is mounted has known dynamics and known external disturbances to its motion, then a dynamic model may be created and updated using, for example, a Kalman filter.

The image sensor need not be one for visible light, but could, for example, be a sensor of infrared (IR), ultraviolet (UV) or other electromagnetic radiation.

The imaging system need not be one that takes discrete snap-shots or “frames”. It could, for example, have detectors that operate continuously in time, such as photo diodes. Time derivatives can be determined from such continuous signals using electrical circuits known in the art.

The imaging system need not produce an actual image output, but instead may produce only the derivatives of brightness in the image needed in the computation of the time to contact, since the image brightness itself does not enter directly into the computation.

The image sensor need not have picture cells on a rectangular grid, but could, for example, be based on a hexagonal or a triangular grid, or even be continuous in space.

This may be of importance for some image sensors that are only responsive to changes in radiation, not being able to measure the absolute irradiance itself.

The calculation of brightness derivatives and sums of products of derivatives and image coordinates need not be performed digitally, but could instead be done in the analog domain, or in a hybrid fashion. Multiplications, for example, could be accomplished using four-quadrant analog multipliers or instead by using multiplying digitial-to-analog (D/A) converters (McQuirk et al 1998).

The summation of sums of products of derivatives and coordinates could be accomplished in the analog domain also, by, for example, summing currents flowing through a resistor or summing charges deposited into a capacitor (McQuirk et al 1998). Analog and hybrid methods for performing the computation may have advantages in terms of component count or fabrication costs.

Claims

1. A time-to-contact estimation device for estimating a time to contact of a moving body to an object in surroundings of the moving body, comprising: an image sensor provided on the moving body, and continuously capturing time-varying images of the surroundings of the moving body;a brightness change calculation unit calculating brightness changes in pixels between the time-varying images captured one after another; anda time-to-contact calculation unit calculating optimized brightness gradients based on the brightness changes calculated by the brightness change calculation unit, and calculating the time to contact of the moving body to the object based on the optimized brightness gradients.

2. The time-to-contact estimation device according to claim 1, wherein the time-to-contact calculation unit calculates the optimized brightness gradients by obtaining a solution that minimizes a sum of the brightness changes in the pixels of at least a portion of the time-varying images.

3. The time-to-contact estimation device according to claim 1, wherein the time-to-contact calculation unit calculates the time to contact of the moving body to the object in different ways depending on whether at least either one of degree of freedom of the object and degree of freedom of the moving body is constrained or none of them are constrained.

4. The time-to-contact estimation device according to claim 1, wherein, when none of degrees of freedom of the object and the moving body are constrained, the time-to-contact calculation unit calculates the time to contact of the moving body to the object by executing an iteration process including plural calculations in which either one of an orientation of the object and a translation path of the moving body is assumed to be fixed in each of the calculations.

5. The time-to-contact estimation device according to claim 3, wherein the time-to-contact calculation unit determines whether at least either one of degree of freedom of the object and degree of freedom of the moving body is constrained or none of them are constrained depending on whether the object is a stationary object or a moving object.

6. The time-to-contact estimation device according to claim 1, wherein the time-to-contact calculation unit calculates the time to contact of the moving body to the object in different ways by classifying cases, depending on an orientation of a planar surface of the object and a translation path along which the moving body approaches the planar surface, into four cases including: a first case in which the planar surface is oriented perpendicular to an optical axis of the image sensor and the moving body approaches the planar surface along the optical axis;a second case in which the planar surface is oriented perpendicular to the optical axis and the moving body approaches the planar surface in an arbitrary direction;a third case in which the planar surface is arbitrarily oriented and the moving body approaches the planar surface along the optical axis; anda fourth case in which the planar surface is arbitrarily oriented and the moving body approaches the planar surface in an arbitrary direction.

7. A method for estimating a time to contact of a moving body to an object in surroundings of the moving body, comprising: an image capturing step of continuously capturing, from the moving body, time-varying images of the surroundings of the moving body;a brightness change calculating step of calculating brightness changes in pixels between the time-varying images captured one after another; anda time-to-contact calculating step of calculating optimized brightness gradients based on the brightness changes calculated in the brightness change calculating step, and calculating the time to contact of the moving body to the object based on the optimized brightness gradients.

8. The method according to claim 7, wherein the time-to-contact calculating step includes, for calculating the optimized brightness gradients, obtaining a solution that minimizes a sum of the brightness changes in the pixels of at least a portion of the time-varying images.

9. The method according to claim 7, wherein, in the time-to-contact calculating step, the time to contact of the moving body to the object is calculated in different ways depending on whether at least either one of degree of freedom of the object and degree of freedom of the moving body is constrained or none of them are constrained.

10. The method according to claim 7, wherein the time-to-contact calculating step includes, when none of degrees of freedom of the object and the moving body are constrained, executing an iteration process including plural calculations in which either one of an orientation of the object and a translation path of the moving body is assumed to be fixed in each of the calculations.

11. The method according to claim 9, wherein, in the time-to-contact calculating step, whether at least either one of degree of freedom of the object and degree of freedom of the moving body is constrained or none of them are constrained is determined depending on whether the object is a stationary object or a moving object.

12. The method according to claim 7, wherein, in the time-to-contact calculating step, the time to contact of the moving body to the object is calculated in different ways by classifying cases, depending on an orientation of a planar surface of the object and a translation path along which the moving body approaches the planar surface, into four cases including: a first case in which the planar surface is oriented perpendicular to an optical axis of the image sensor and the moving body approaches the planar surface along the optical axis;a second case in which the planar surface is oriented perpendicular to the optical axis and the moving body approaches the planar surface in an arbitrary direction;a third case in which the planar surface is arbitrarily oriented and the moving body approaches the planar surface along the optical axis; anda fourth case in which the planar surface is arbitrarily oriented and the moving body approaches the planar surface in an arbitrary direction.