SYSTEMS AND METHODS FOR ESTIMATING THE STRUCTURE AND MOTION OF AN OBJECT

[0001] This application claims priority to co-pending U.S. Provisional Application Ser. No. 61/496,191, filed Jun. 13, 2011, which is hereby incorporated by reference herein in its entirety.

[0002] This invention was made with government support under the National Science Foundation (NSF) Career Award Number 0547448, NSF Award Number 0901491, and the Department of Energy Grant Number DE-FG04-86NE37967. The government has certain rights in the invention.

[0003] In certain situations, it is useful to be able to estimate the structure and motion of an object. For example, it is desirable to estimate the structure and relative motion of stationary or moving objects that an unmanned aerial vehicle (UAV) must avoid. In prior solutions to this problem, knowledge of the Euclidean geometry of the objects or full knowledge of the camera motion of the UAV is required. It would be desirable to have a system and method that do not require that knowledge.

[0015] As described above, it would be desirable to be able to estimate the structure and motion of objects without knowledge of the Euclidean geometry of the objects or full knowledge of the camera motion. Disclosed herein are systems and methods in which such structure and motion can be estimated. In some embodiments, the systems and methods can be used to determine the structure and motion of stationary objects. In other embodiments, the systems and methods can be used to determine the structure and motion of stationary and/or moving objects.

[0016] In a first embodiment, a reduced order nonlinear observer is proposed for the problem of structure and motion (SaM) estimation of a stationary object observed by a moving calibrated camera. In comparison to existing work which requires some knowledge of the Euclidean geometry of an observed object or full knowledge of the camera motion, the developed reduced order observer only requires one camera linear velocity and corresponding acceleration to asymptotically identify the Euclidean coordinates of the feature points attached to an object (with proper scale reconstruction) and the remaining camera velocities. The unknown linear velocities are assumed to be generated using a model with unknown parameters. The unknown angular velocities are determined from a robust estimator which uses a standard homography decomposition algorithm applied to tracked feature points. A Lyapunov analysis is provided to prove the observer asymptotically estimates the unknown states under a persistency of excitation (PE) condition.

[0017] An objective of the classic “structure from motion (SfM)” problem is to estimate the Euclidean coordinates of tracked feature points attached to an object (i.e., three-dimensional (3D) structure) provided the relative motion between the camera and the object is known. The converse of the SfM problem is the “motion from structure (MfS)” problem where the relative motion between the camera and the object is estimated based on known geometry of the tracked feature points attached to an object. An extended problem disclosed herein is a “structure and motion (SaM)” problem where the objective is to estimate the Euclidean geometry of the tracked feature points as well as the relative motion between the camera and tracked feature points. The SaM problem is a fundamental problem and some examples indicate that SaM estimation is only possible up to a scale when a pinhole camera model is used. To recover the scale information, either the linear camera velocities or partial information about the structure of the object, e.g. a known length between two feature points in a scene, is required. In this disclosure, a partial solution to the SaM problem is presented, where the objective is to asymptotically identify the Euclidean 3D coordinates of tracked feature points and the camera motion, provided at least one linear velocity of the camera is known. The angular velocity is estimated using a filter. The estimated angular velocity and a measured linear velocity are combined to estimate the scaled 3D coordinates of the feature points.

[0018] Solutions to the SfM, MfS, and the SaM problems can be broadly classified as offline methods (batch methods) and online methods (iterative methods). Batch methods extract an image data set from a given image sequence and then the 3D structure is reconstructed from the data set. These methods are usually based on nonlinear optimization, projective methods, or invariant-based methods. Often, batch methods lack an analytical analysis of convergence. The main drawback of batch methods is that they cannot be used to execute online/real-time tasks. Thus, the need arises for iterative or online methods with analytical guarantees of convergence.

[0019] Online methods typically formulate the SfM and MfS problems as a continuous differential equation, where the image dynamics are derived from a continuous image sequence. Online methods often rely on the use of an Extended Kalman filter (EKF). Kalman filter based approaches also lack a convergence guarantee and can converge to wrong solutions in practical scenarios. Also, a priori knowledge about the noise is required for such solutions. In comparison to Kalman filter-based approaches, some researchers have developed nonlinear observers for SfM with analytical proofs of stability. For example, a discontinuous high-gain observer called identifier-based observer (IBO) has been evaluated under the assumption of known camera motion. In one solution, a discontinuous sliding-mode observer is developed which guarantees exponential convergence of the states to an arbitrarily small neighborhood, i.e., uniformly ultimately bounded (UUB) result. A continuous observer which guarantees asymptotic structure estimation has been evaluated under the assumption of known camera motion. An asymptotically stable reduced-order observer has also been evaluated to estimate the structure given known camera motion. Under the assumption that a known Euclidean distance between two feature points is known, a nonlinear observer can be used to asymptotically identify the camera motion. In contrast to these approaches, the methods described in this disclosure do not assume any model knowledge of the structure and only require one linear velocity and the corresponding acceleration.

[0020] Various batch and iterative methods have been developed to solve the SaM problem upto a scale. However, in comparison to SfM and MfS results, sparse literature is available where the SaM problem is formulated in terms of continuous image dynamics with associated analytical stability analysis. In one study, an algorithm was presented to estimate the structure and motion parameters up to a scaling factor. In another study, a perspective realization theory for the estimation of the shape and motion of a moving planar object observed using a static camera up to a scale was evaluated. Recently, a nonlinear observer was developed to asymptotically identify the structure given the camera motion (i.e., the SfM problem) or to asymptotically identify the structure and the unknown time-varying angular velocities given all three linear velocities. In a further study structure and linear velocities were estimated given partial structure information such as length between two points on an object, which may be difficult in practice for random objects. In another recent study, the IBO or IBO approach was used to estimate the structure and the constant angular velocity of the camera given all three linear velocities, which not be possible in practical scenarios such as a camera attached to a unmanned vehicle where sideslip velocities may not be available. The problem of estimating structure, time varying angular velocities, and time varying linear velocities of the camera without knowledge of partial structure information remains an unsolved problem.

[0021] The technical challenge presented by the SaM problem is that the image dynamics are scaled by an unknown factor, and the unknown structure is multiplied by unknown motion parameters. As described in the following discussions, a challenge is to estimate a state in the open loop dynamics that appears nonlinearly inside a matrix that is multiplied by a vector of unknown linear and angular velocity terms. By assuming that the velocities are known, or some model knowledge exists, previous online efforts have been able to avoid the problem of separately estimating multiplicative uncertainties. One contribution of this disclosure is a strategy to segregate the multiplicative uncertainties, and then to develop a reduced order nonlinear observer to address the SaM problem where the structure (i.e., the properly scaled relative Euclidean coordinates of tracked feature points), the time-varying angular velocities, and two unknown time-varying linear velocities are estimated (i.e., one relative linear velocity is assumed to be known along with a corresponding acceleration). The result exploits an uncertain locally Lipschitz model of the unknown linear velocities of the camera. The strategic use of a standard homography decomposition is used to estimate the angular velocities, provided the intrinsic camera calibration parameters are known and feature points can be tracked between images. A persistency of excitation (PE) condition is formulated, which provides an observability condition that can be physically interpreted as the known camera linear velocity should not be zero over any small interval of time, and the camera should not be moving along the projected ray of a point being tracked. A Lyapunov-based analysis is provided that indicates the SaM observer errors are globally asymptotically regulated provided the PE condition is satisfied. By developing a reduced order observer to segregate and estimate the multiplicative uncertainties, new applications can be addressed including: range and velocity estimation using a camera fixed to a moving vehicle where only the forward velocity/acceleration of the vehicle is known, range and velocity estimation using an unmanned air vehicle (UAV) using only a forward velocity/acceleration sensors, etc.

[0022] Consider a moving camera that views four or more planar and non-collinear feature points (denoted by j={1, 2, . . . , n} ∀n≧4) lying fixed in a visible plane π_r, attached to an object in front of the camera. Let F_r be a static coordinate frame attached to the object. A static reference orthogonal coordinate frame F_c* is attached to the camera at the location corresponding to an initial point in time t₀ where the object is in the camera field of view (FOV). After the initial time, an orthogonal coordinate frame F_c attached to the camera undergoes some rotation R(t)εSO(3) and translation x_ƒ(t)ε³ away from FF_c*.

[0023] The Euclidean coordinates m_j(t)ε³ of the feature points expressed in the current camera frame F_c and the respective normalized Euclidean coordinates m_j(t)ε³ are defined as

$\begin{array}{cc}{\stackrel{\_}{m}}_j\left(t\right)={\left[x_{1j}\left(t\right)x_{2j}\left(t\right)x_{3j}\left(t\right)\right]}^T,m_j\left(t\right)={\left[\frac{x_{1j}\left(t\right)}{x_{3j}\left(t\right)}\frac{x_{2j}\left(t\right)}{x_{3j}\left(t\right)}1\right]}^T.& \left(1\right)\end{array}$

The constant Euclidean coordinates and the normalized coordinates of the feature points expressed in the camera frame F_c* are denoted by m_j*ε³, and m_j*ε³ respectively and are given by Equation (1) superscripted by a ‘*’. Consider a closed and bounded set y⊂³. Auxiliary state vectors y_j*=[y_1j*, y_2j*, y_3j*]^Tεy and y_j(t)=[y_1j(t), y_2j(t), y_3j]^Tεy are constructed from Equation (1) as

$\begin{array}{cc}{y}_j^{*}={\left[\frac{x_{1j}^{*}}{x_{3j}^{*}}\frac{x_{2j}^{*}}{x_{3j}^{*}}\frac{1}{x_{3j}^{*}}\right]}^T,y_j={\left[\frac{x_{1j}}{x_{3j}}\frac{x_{2j}}{x_{3j}}\frac{1}{x_{3j}}\right]}^T.& \left(2\right)\end{array}$

The corresponding feature points m_j* and m_j(t) viewed by the camera from two different locations (and two different instances in time) are related by a depth ratio

${\alpha }_j\left(t\right)\stackrel{\Delta }{=}\frac{x_{3j}^{*}}{x_{3j}}\in ℝ$

and a homography matrix H(t)ε^3×3 as

m_j=α_jHm_j*.  (3)

The homography matrix can be decomposed to obtain a rotation and scaled translation of the camera between two views. The decomposition of the homography leads to two solutions, one of which is physically relevant. Using projective geometry, the normalized Euclidean coordinates m_j* and m_j(t) can be related to the pixel coordinates in the image space as

p_j=Am_j,p_j*=Am_j*  (4)

where p_j(t)=[u_j, ν_j, 1]^T is a vector of the image-space feature point coordinates u_j(t), ν_j(t)ε defined on the closed and bounded set I⊂³, and Aε^3×3 is a constant, known, invertible intrinsic camera calibration matrix. Since A is known, Equation (4) can be used to recover m_j(t), which can be used to partially reconstruct the state y(t).

[0024] Assumption 1:

[0025] The relative Euclidean distance x_3j(t) between the camera and the feature points observed on the target is upper and lower bounded by some known positive constants (i.e., the object remains within some finite distance away from the camera).

[0026] Remark 1:

[0027] The definition in Equation (2) along with Assumption 1 can be used to conclude that y₃ (t) remains in a set Ω defined as Ω={y₃|y₃≦y₃≦ y₃} where y₃, y₃ε denote known positive bounding constants. Likewise, since the image coordinates are constrained (i.e., the target remains in the camera field of view) the relationships in Equations (1), (2), and (4) along with the fact that A is invertible can be used to conclude that y₁≧|y_1j(t)|≧y₁, y≧|y_2j(t)|≧y₂ where y₁, y₂, y₁, y₂ε denote known positive bounding constants. Furthermore, the velocity parameters of the camera are bounded by constants.

[0028] At some spatiotemporal instant, the camera views a point q on the object. The point q can be expressed in the coordinate system F_c as

m=x_ƒ+Rx_o′q  (5)

where x_o′q is a vector from the origin of the static coordinate system F_r, attached to the object at the point q. By differentiating Equation (5), the relative motion of q can be expressed as

m=[ω]_x m+b

where m(t) is defined in Equation (1), [ω]_x, ε_3×3 denotes a skew symmetric matrix formed from the angular velocity vector of the camera ω(t)=[ω₁, ω₂, ω₃]^T εWω₃^T εW, and b(t)=[b₁, b₂, b₃]^TεB denotes the linear velocity of the camera. The sets W and B are closed and bounded sets such that W⊂³ and B⊂³. Using Equations (2) and (5), the dynamics of the partially measurable state y(t) can be expressed as

$\begin{array}{cc}\left[\begin{array}{c}{\stackrel{.}{y}}_1\\ {\stackrel{.}{y}}_2\\ {\stackrel{.}{y}}_3\end{array}\right]=\left[\begin{array}{cccccc}{y}_3& 0& -y_1y_3& -y_1y_2& 1+y_1^2& -y_2\\ 0& y_3& -y_2y_3& -\left(1+y_2^2\right)& y_1y_2& y_1\\ 0& 0& y_3^2& y_2y_3& -y_1y_3& 0\end{array}\right]\left[\begin{array}{c}b\\ \omega \end{array}\right]& \left(6\right)\end{array}$

where y₁(t) and y₂(t) can be measured through the transformation in Equation (4).

[0029] An estimator can be designed for the perspective dynamic system in Equation (6), where the angular velocity is considered unknown and only one of the linear velocities (i.e., b₃) and respective acceleration (i.e., {dot over (b)}₃) is available. Moreover, an uncertain dynamic model of the linear velocity b(t) is assumed to be available as

{dot over (b)}_i(t)=q(b_i,t)b_i,∀i={1,2}  (7)

where q(b_i,t)ε is a known locally Lipschitz function of unknown states. To facilitate the design and analysis of the subsequent observer, a new state

$u\left(t\right)\in U⋐{ℝ}^2\stackrel{\Delta }{=}{\left[u_1\left(t\right)=y_3b_1,u_2\left(t\right)=y_3b_2\right]}^T,$

is defined where U is a closed and bounded set. After utilizing Equations (6) and (7), the dynamics for u₁(t), u₂(t) can be expressed as

{dot over (u)}_i=y₃b₃u_i+(y₂ω₁−y₁ω₂)u_i+q(b_i)u_i,∀i={1,2}.  (8)

From Equations (6) and (8) the dynamics of the known states y₁(t), y₂(t) and the unknown state θ(t)=[y₃ u₁ u₂]^T are

$\begin{array}{cc}\left[\begin{array}{c}{\stackrel{.}{y}}_1\\ {\stackrel{.}{y}}_2\end{array}\right]=\underset{\underset{J\left(y,b_3\right)}{}}{\left[\begin{array}{ccc}-y_1b_3& 1& 0\\ -y_2b_3& 0& 1\end{array}\right]}\left[\begin{array}{c}{y}_3\\ u_1\\ u_2\end{array}\right]+\underset{\underset{\psi \left(y,\omega \right)}{}}{\left[\begin{array}{c}-y_1y_2{\omega }_1+\left(1+y_1^2\right){\omega }_2-y_2{\omega }_3\\ -\left(1+y_2^2\right){\omega }_1+y_1y_2{\omega }_2+y_1{\omega }_3\end{array}\right]}\mathrm{and}& \left(9\right)\\ \left[\begin{array}{c}{\stackrel{.}{y}}_3\\ {\stackrel{.}{u}}_1\\ {\stackrel{.}{u}}_2\end{array}\right]=\underset{\underset{g\left(\theta ,\omega ,y_1,y_2,b_3\right)}{}}{\left[\begin{array}{c}{y}_3^2b_3+\left(y_2{\omega }_1-y_1{\omega }_2\right)y_3\\ y_3b_3u_1+\left(y_2{\omega }_1-y_1{\omega }_2\right)u_1+q\left(b_1\right)u_1\\ y_3b_3u_2+\left(y_2{\omega }_1-y_1{\omega }_2\right)u_2+q\left(b_2\right)u_2\end{array}\right]}.& \left(10\right)\end{array}$

Since y₁(t) and y₂(t) are measurable, from Equations (1) and (4) the Euclidean structure m(t) can be estimated once the state y₃(t) is determined. Since the dynamics of the outputs y₁(t), y₂ (t) are affine in the unknown state θ(t), a reduced order observer can be developed based on this relationship for the unknown state θ(t). The subsequent development is based on the strategy of constructing the estimates

$\hat{\theta }\left(t\right)\stackrel{\Delta }{=}{\left[{\hat{y}}_3{\hat{u}}_1{\hat{u}}_2\right]}^T\in {ℝ}^3.$

To quantify the SaM estimation objective, an estimation error {tilde over (θ)}(t)=[{tilde over (θ)}₁ {tilde over (θ)}₂ {tilde over (θ)}₃]^Tε³ is defined as

$\begin{array}{cc}\tilde{\theta }\left(t\right)\stackrel{\Delta }{=}{\left[y_3-{\hat{y}}_3u_1-{\hat{u}}_1u_2-{\hat{u}}_2\right]}^T.& \left(11\right)\end{array}$

[0030] Assumption 2:

[0031] The function q(b_i, t) ∀i={1,2} is locally Lipschitz where q(b₁)−q({circumflex over (b)}₁)=λ₁(b₁−{circumflex over (b)}₁) and q(b₂)−q({circumflex over (b)}₂)=λ₂(b₂−{circumflex over (b)}₂) where λ₁ and λ₂ are Lipschitz constants.

[0032] Remark 2:

[0033] The linear velocity model in Equation (7) is restricted to motions that are satisfied by Assumption 2; yet, various classes of trajectories satisfy this assumption (e.g., straight line trajectories, circles, some periodic trajectories, etc.).

[0034] Assumption 3:

[0035] The function J(y₁, y₂, b₃) defined in Equation (9) satisfies the persistency of excitation condition, i.e., ∃γ, δ>0:∫_t^t+δJ^T(y₁(τ), y₂(τ))J(y₁(τ), b₃(τ)dτ≧γI∀t≧0.

[0036] Remark 3:

[0037] Assumption 3 is violated if ∃t₁|∃t>t₁, b₃ (t)=0 or y₁(t)=c₁, y₂(t)=c₂. That is, Assumption 3 is valid unless there exists a time t₁ such that for all t>t₁ the camera translates along the projected ray of an observed feature point.

[0038] Assumption 4:

[0039] The linear camera velocities b(t) are upper and lower bounded by constants.

[0040] Remark 4:

[0041] The following bounds can be developed using Assumption 1, Remark 1 and the definitions of u₁(t) and u₂ (t)

u_1min≦u₁≦u_1max,u_2min≦u₂≦u_2max.

[0042] Step I:

[0043] Angular Velocity Estimation Solutions are known that can be used to determine the relative angular velocity between the camera and a target. To quantify the rotation mismatch between F_c* and F_c, a rotation error vector e_ω(t)ε is defined by the angle-axis representation as

$\begin{array}{cc}{e}_{\omega }\stackrel{\Delta }{=}u_{\omega }\left(t\right){\theta }_{\omega }\left(t\right)& \left(12\right)\end{array}$

where u_ω(t)ε³ represents a unit rotation axis, and θ_ω(t)ε denotes the rotation angle about u_ω(t) that is assumed to be confined to region −π<θ_ω(t)<π. The angle θ_ω(t) and axis u_ω(t) can be computed using the rotation matrix R(t) obtained by decomposing the homography matrix H(t) given by the relation in Equation (3). Taking time derivative of Equation (12) yields

ė_ω=L_ωω  (13)

where L_ω(t)εR^3×3 denotes an invertible Jacobian matrix. A robust integral of the sign of the error (RISE)-based observer ê_ω(t)ε³ is generated as

ê_ω=(K_ω+I_3×3){tilde over (e)}_ω(t)∫_t0^t(K_ω+I_3×3){tilde over (e)}_ωdτ+ν

{dot over (ν)}=ρ_ωsgn({tilde over (e)}_ω)  (14)

where K_ω, ρ_ωε^3×3 positive constant diagonal gain matrices, and {tilde over (e)}_ω(t)ε³ quantifies the observer error as

${\tilde{e}}_{\omega }\left(t\right)\stackrel{\Delta }{=}e_{\omega }-{\hat{e}}_{\omega }.$

A Lyapunov-based stability analysis proves

{circumflex over (ė)}_ω(t)−ė_ω(t)→0 as t→∞  (15)

and that all closed-loop signals are bounded. Based on Equations (13) and (15), the angular velocity can be determined as

{circumflex over (ω)}(t)=L_ω⁻¹{circumflex over (ė)}_ω(t)as t→∞  (16)

An angular velocity estimation error

$\tilde{\omega }\left(t\right)\in {ℝ}^3\stackrel{\Delta }{=}[{\tilde{\omega }}_1\left(t\right),$

{tilde over (ω)}₂(t), {tilde over (ω)}₃(t)]^T is defined as ω_i(t)=ω_i (t)−{circumflex over (ω)}_i(t), ∀i={1, 2, 3}. The angular velocity estimator given by Equation (14) is asymptotically stable; thus, the angular velocity estimation error ∥ ω(t)∥→0 as t→∞.

[0044] Step II:

[0045] Structure Estimation: A reduced order observer for θ(t) is designed as

$\begin{array}{cc}\left[\begin{array}{c}{\hat{y}}_3\\ {\hat{u}}_1\\ {\hat{u}}_2\end{array}\right]=\left[\begin{array}{c}{\stackrel{\_}{y}}_3\\ {\stackrel{\_}{u}}_1\\ {\stackrel{\_}{u}}_2\end{array}\right]+\Gamma \left[\begin{array}{c}\frac{-b_3\left(y_1^2+y_2^2\right)}{2}\\ y_1\\ y_2\end{array}\right]& \left(17\right)\end{array}$

where the state vector [ y₃ ū₁ ū₂]^T is updated using the following update law

$\begin{array}{cc}\left[\begin{array}{c}{\stackrel{.}{\stackrel{\_}{y}}}_3\\ {\stackrel{.}{\stackrel{\_}{u}}}_1\\ {\stackrel{\_}{\stackrel{.}{u}}}_2\end{array}\right]=\underset{\underset{g\left(\hat{\theta },\hat{\omega },y_1,y_2,b_3\right)}{}}{\left[\begin{array}{c}{\hat{y}}_3^2b_3+\left(y_2{\hat{\omega }}_1-y_1{\hat{\omega }}_2\right){\hat{y}}_3\\ {\hat{y}}_3b_3{\hat{u}}_1+\left(y_2{\hat{\omega }}_1-y_1{\hat{\omega }}_2\right){\hat{u}}_1+q\left({\hat{b}}_1\right){\hat{u}}_1\\ {\hat{y}}_3b_3{\hat{u}}_2+\left(y_2{\hat{\omega }}_1-y_1{\hat{\omega }}_2\right){\hat{u}}_2+q\left({\hat{b}}_2\right){\hat{u}}_2\end{array}\right]}-\Gamma J^T\left(\left[\begin{array}{ccc}-y_1b_3& 1& 0\\ -y_2b_3& 0& 1\end{array}\right]\left[\begin{array}{c}{\hat{y}}_3\\ {\hat{u}}_1\\ {\hat{u}}_2\end{array}\right]+\underset{\underset{\psi \left(y,\tilde{\omega }\right)}{}}{\left[\begin{array}{c}-y_1y_2{\hat{\omega }}_1+\left(1+y_1^2\right){\hat{\omega }}_2-y_2{\hat{\omega }}_3\\ -\left(1+y_2^2\right){\hat{\omega }}_1+y_1y_2{\hat{\omega }}_2+y_1{\hat{\omega }}_3\end{array}\right]}\right)+\Gamma \left[\begin{array}{c}\frac{b_3\left(y_1^2+y_2^2\right)}{2}\\ 0\\ 0\end{array}\right].& \left(18\right)\end{array}$

In Equation (18), Γε⁺, {circumflex over (ω)}_i(t) are given by Equation (16), and J(y,b₃) is defined in Equation (9). Differentiating Equation (11) and using Equations (9), (10), (17) and (18) yields the following closed-loop observer error dynamics

$\left[\begin{array}{c}{\tilde{\theta }}_1\\ {\tilde{\theta }}_2\\ {\tilde{\theta }}_3\end{array}\right]=\left[\begin{array}{c}\left(y_3+{\hat{y}}_3\right)b_3{\tilde{\theta }}_1+\left(y_2{\omega }_1-y_1{\omega }_2\right){\tilde{\theta }}_1+\left(y_2{\tilde{\omega }}_1-y_1{\tilde{\omega }}_2\right){\hat{y}}_3\\ y_3b_3{\tilde{\theta }}_2+b_3{\hat{u}}_1{\tilde{\theta }}_1+\left(y_2{\hat{\omega }}_1-y_1{\hat{\omega }}_2\right){\tilde{\theta }}_2+\left(y_2{\tilde{\omega }}_1-y_1{\tilde{\omega }}_2\right){\hat{u}}_1+q\left(b_1\right){\tilde{\theta }}_2+{\lambda }_1\left(b_1-{\hat{b}}_1\right){\hat{u}}_1\\ y_3b_3{\tilde{\theta }}_3+b_3{\hat{u}}_2{\tilde{\theta }}_1+\left(y_2{\hat{\omega }}_1-y_1{\hat{\omega }}_2\right){\tilde{\theta }}_3+\left(y_2{\tilde{\omega }}_1-y_1{\tilde{\omega }}_2\right){\hat{u}}_2+q\left(b_2\right){\tilde{\theta }}_3+{\lambda }_2\left(b_2-{\hat{b}}_2\right){\hat{u}}_2\end{array}\right]+\Gamma J^T\left(\left[\begin{array}{ccc}-y_1b_3& 1& 0\\ -y_2b_3& 0& 1\end{array}\right]\left[\begin{array}{c}{\hat{y}}_3\\ {\hat{u}}_1\\ {\hat{u}}_2\end{array}\right]+\left[\begin{array}{c}-y_1y_2{\hat{\omega }}_1+\left(1+y_1^2\right){\hat{\omega }}_2-y_2{\hat{\omega }}_3\\ -\left(1+y_2^2\right){\hat{\omega }}_1+y_1y_2{\hat{\omega }}_2+y_1{\hat{\omega }}_3\end{array}\right]-\left[\begin{array}{c}{\stackrel{.}{y}}_1\\ {\stackrel{.}{y}}_2\end{array}\right]\right).$

Using the output dynamics from Equation (9), the error dynamics can be rewritten as

θ=g(θ,ωy₁,y₂,b₃)−g({circumflex over (θ)},{circumflex over (ω)},y₁,y₂,b₃)−ΓJ^T(J{tilde over (θ)}+ψ(y, ω)).

Using Assumption 2, the following relationship can be developed

g(θ,ω,y₁,y₂,b₃)−g({circumflex over (θ)}{circumflex over (ω)}y₁,y₂b₃)=Λ θ+α(y₁,y₂,{circumflex over (θ)}{tilde over (ω)})  (19)

where

$\begin{array}{cc}\Lambda \left(t\right)=\left[\begin{array}{ccc}\left(y_3+{\hat{y}}_3\right)b_3+\left(y_2{\omega }_1-y_1{\omega }_2\right)& 0& 0\\ b_3{\hat{u}}_1+\frac{{\lambda }_1{\hat{u}}_1^2}{y_3{\hat{y}}_3}& y_3b_3+q\left(b_1\right)+\left(y_2{\hat{\omega }}_1-y_1{\hat{\omega }}_2\right)+\frac{{\lambda }_1{\hat{u}}_1}{y_3}& 0\\ b_3{\hat{u}}_2+\frac{{\lambda }_2{\hat{u}}_2^2}{y_3{\hat{y}}_3}& 0& y_3b_3+q\left(b_2\right)+\left(y_2{\hat{\omega }}_1-y_1{\hat{\omega }}_2\right)+\frac{{\lambda }_2{\hat{u}}_2}{y_3}\end{array}\right],& \left(20\right)\end{array}$

such that sup_t∥Λ(t)∥≦ρ, and

$\begin{array}{cc}\alpha \left(y_1,y_2,\hat{\theta },\tilde{\omega }\right)=\left[\begin{array}{c}\left(y_2{\tilde{\omega }}_1-y_1{\tilde{\omega }}_2\right){\hat{y}}_3\\ \left(y_2{\tilde{\omega }}_1-y_1{\tilde{\omega }}_2\right){\hat{u}}_1\\ \left(y_2{\tilde{\omega }}_1-y_1{\tilde{\omega }}_2\right){\hat{u}}_2\end{array}\right].& \left(21\right)\end{array}$

Using Equation (19) the error dynamics can be written as

θ=Λ θ−ΓJ^TJ θ−ΓJ^TΨ(y, ω)+α(y₁,y₂{circumflex over (θ)} ω).  (22)

[0046] The results from the angular velocity estimator in Section IV-A1 prove that ω₁, ω₂, ω₃→0 therefore, Equation (21) can be used to conclude that Ψ(•)→0 and α(•)→0 as t→∞.

[0047] Theorem 1:

[0048] If Assumptions 1-4 are satisfied, the reduced order observer in Equations (17) and (18) asymptotically estimates θ(t) in the sense that ∥ θ(t)∥→0 as t→∞.

[0049] Proof:

[0050] The stability of the error system in Equation (22) can be proved using a converse Lyapunov theorem. Consider the nominal system

θ=ƒ( θ)=−ΓJ^TJ θ.  (23)

The error system in Equation (23) is globally exponentially stable if Assumption 3 is satisfied. Hence, {tilde over (θ)}(t) satisfies the inequality ∥{tilde over (θ)}(t)∥≦∥{tilde over (θ)}(t₀)∥α₁ e^{−α2(t−t0)}, where α₁, α₂ε⁺, and α₂ is directly proportional to Γ and inversely proportional to δ. Consider a set D={ θ(t)ε³|∥{tilde over (θ)}(t)∥<∞}. Using a converse Lyapunov theorem there exists a function V:[0,∞)×D→ that satisfies

$\begin{array}{cc}{c}_1{\tilde{\theta }\left(t\right)}^2\le V\left(t,\tilde{\theta }\right)\le c_2{\tilde{\theta }\left(t\right)}^2,\frac{\partial V}{\partial t}+\frac{\partial V}{\partial \tilde{\theta }}\left(-\Gamma J^TJ\tilde{\theta }\right)\le -c_3{\tilde{\theta }\left(t\right)}^2,\frac{\partial V}{\partial \tilde{\theta }}\le c_4\tilde{\theta }\left(t\right)& \left(24\right)\end{array}$

for some positive constants c₁, c₂, c₃, c₄. Using V(t, {tilde over (θ)}) as a Lyapunov function candidate for the perturbed system in Equation (22), the derivative of V(t,{tilde over (θ)}) along the trajectories of Equation (22) is given by

$\stackrel{.}{V}\left(t,\tilde{\theta }\right)=\frac{\partial V}{\partial t}+\frac{\partial V}{\partial \tilde{\theta }}\left(-\Gamma J^TJ\tilde{\theta }\right)+\frac{\partial V}{\partial \tilde{\theta }}\left(\left(\Lambda \right)\tilde{\theta }\right)+\frac{\partial V}{\partial \tilde{\theta }}\left(-\Gamma J^T\Psi \left(y,\tilde{\omega }\right)+\alpha \left(y_1,y_2,\hat{\theta },\tilde{\omega }\right)\right).$

Using the bounds in Equation (24) the following inequality is developed

{dot over (V)}(t,{tilde over (θ)})≦−(c₃−c₄ρ)∥{tilde over (θ)}∥²+c₄d∥{tilde over (θ)}∥  (25)

where ρ is an upper bound on (20) and d(t)=∥Γ∥∥J^T∥∥Ψ(y,{tilde over (ω)})∥+∥α(y₁,y₂,{circumflex over (θ)}{tilde over (ω)})∥ where d(t)→0 as t→∞. The estimates of c₃ and c₄ are given by

$c_3=\frac{1}{2},c_4=\frac{2{\alpha }_1}{\left({\alpha }_2-L\right)}\left[1-{}^{-\frac{\left({\alpha }_2-L\right)\ln \left(2{\alpha }_1^2\right)}{2{\alpha }_2}}\right]$

where Lε⁺ is an upper bound on the norm of Jacobian matrix

$\frac{\partial f\left(\tilde{\theta }\right)}{\partial \tilde{\theta }},$

where ƒ({tilde over (θ)}) is defined in Equation (23). Since α₂ is directly proportional to the gain Γ, the inequality

$c_3-c_4\rho =\frac{1}{2}-\frac{2{\alpha }_1\rho }{\left({\alpha }_2-L\right)}\left[1-{}^{-\frac{\left({\alpha }_2-L\right)\ln \left(2{\alpha }_1^2\right)}{2{\alpha }_2}}\right]>0$

can be achieved by choosing the gain Γ sufficiently large. Using Equations (24) and (25), the following bound is obtained

$\begin{array}{cc}\tilde{\theta }\left(t\right)\le \sqrt{\frac{c_2}{c_1}}{}^{\frac{c_4}{2c_1}}\tilde{\theta }\left(t_0\right){}^{-\beta \left(t-t_0\right)}+\frac{c_4}{2c_1}{}^{\frac{c_4}{2c_1}}{\int }_0^t{}^{-\beta \left(t-\tau \right)}\left(\tau \right)\tau & \left(26\right)\end{array}$

where a constant convergence rate β>0 can be increased by increasing c₃. From Equation (26), ∥{tilde over (θ)}(t)∥εL_∞, thus {tilde over (θ)}₁(t), {tilde over (θ)}₂ (t), {tilde over (θ)}₃(t)εL_∞. Since {tilde over (θ)}₁(t), {tilde over (θ)}₂ (t), {tilde over (θ)}₃(t)εL_∞, and the fact that θ₁(t), θ₂(t), θ₃(t)εL_∞ can be used to conclude that {circumflex over (θ)}₁(t), {circumflex over (θ)}₂(t), {circumflex over (θ)}₃(t)εL_∞. Using the result from Section IV-A1 that ∥{tilde over (ω)}(t)∥→0 as t→∞, the functions ∥Ψ(y,{tilde over (ω)})∥, ∥α(y₁, y₂,{circumflex over (θ)}, {tilde over (ω)})∥→0 as t→∞. Hence, d(t)→0 as t→∞ and d(t)εLL_∞. Since d(t)→0 as t→∞ and d(t)εLL_∞, by the Lebesgue dominated convergence theorem

$\begin{array}{c}{\lim }_{t->\infty }{\int }_0^t{}^{-\beta \left(t-\tau \right)}\left(\tau \right)\tau ={\int }_0^t{}^{-\mathrm{\beta \sigma }}{\lim }_{t->\infty }\left(\tau -\sigma \right)\sigma \\ =\frac{{\lim }_{t->\infty }d\left(t\right)}{\beta }\\ =0.\end{array}$

It can be shown that ∥{tilde over (θ)}(t)∥→0 as t→∞. Hence, the reduced order estimator in Equations (17) and (18) identifies the structure of observed feature points and unknown camera motion asymptotically. Since y₃(t), u₁(t), and u₂(t) can be estimated, the motion parameters b₁(t) and b₂(t) can be recovered based on the definition of u(t).

[0051] As described above, a reduced order observer is developed for the estimation of the structure (i.e. range to the target and Euclidean coordinates of the feature points) of a stationary target with respect to a moving camera, along with two unknown time-varying linear velocities and the angular velocity. The angular velocity is estimated using homography relationships between two camera views. The observer requires the image coordinates of the points, a single linear camera velocity, and the corresponding linear camera acceleration in any one of the three camera coordinate axes. Under a physically motivated PE condition, asymptotic convergence of the observer is guaranteed. Having at least some partial motion knowledge and a similar observability condition to the given PE condition are likely necessary in future solutions to the SaM problem. However, future efforts could potentially eliminate the need for any model of the vehicle trajectory (even if uncertain as in this result) and eliminate the need for an acceleration measurement.

[0052] FIGS. 1A-1C illustrate an example method for estimating the structure and motion of a stationary object consistent with the above discussion. Beginning with block 100, the rotation error vector ê_ω and the state vector [ y₃, ū₁, ū₂] are initialized. As can be appreciated from the foregoing discussion, the rotation error vector ê_ω is a representation of the rotation between a previous (initial) image or frame and a current image or frame, y₃ is an auxiliary signal for the inverse depth of the feature point (ŷ₃), ū₁, is an estimate of the camera linear velocity in the x direction multiplied by the inverse of the depth, and ū₂ is an estimate of the camera linear velocity in the y direction multiplied by the inverse of the depth. Substantially any initial values (e.g., random numbers) can be used to initialize those vectors, although the closer the selected values are to the actual values that are being estimated, the faster convergence will be attained. In some embodiments, the values that are selected depend upon the underlying conditions in which the estimation is being performed.

[0053] Once the rotation error vector and the state vector have been initialized, an initial image of an object can be captured using a moving camera, as indicated in block 102. FIG. 2 schematically illustrates such image capture with a camera. By way of example, the object is a stationary terrestrial object and the camera is mounted to a moving object, such as an unmanned aerial vehicle (UAV). Once the initial image is captured, feature points of the object within the image are extracted, as indicated in block 104. By way of example, four coplanar feature points are extracted or eight non-coplanar feature points are extracted. In some embodiments, coplanar points can be identified using image segmentation techniques. With reference to block 106, the x and y pixel coordinates of each image feature point are recorded and, with reference to block 108, the linear velocity and linear acceleration of the camera motion at the time of image capture, which are measured by sensors (e.g., accelerometers) associated with the camera, are recorded. In some embodiments, the linear velocity and acceleration along the z direction are recorded. In alternative embodiments, however, the linear velocity and acceleration along the x or y direction can be recorded. At this point, the Euclidean coordinates m_j of the feature points are computed and m_j is designated as m_j*, as indicated in block 110. As can be appreciated from the foregoing discussion, m_j is a transformed image point in the camera frame that is defined by Equation (4). That equation indicates that m, can be determined from a pixel location p_j and a camera calibration matrix A.

[0054] Referring next to block 112, the next image of the object is captured using the moving camera and, as before, feature points within the image are extracted (block 114), the x and y pixel coordinates of each image feature point is recorded (block 116 of FIG. 1B), and the linear velocity and linear acceleration of the camera motion at the time of image capture is recorded (block 118). At this point, the Euclidean coordinates m_j of the feature points are computed using the image pixel coordinates and Equation (4), as indicated in block 120.

[0055] Turning to block 122, the homography matrix H is estimated using Equation (3) and it is decomposed to estimate the rotation matrix R(t) between the initial camera coordinate frame and the current camera coordinate frame. As can be appreciated from the foregoing discussion, the homography matrix H is a geometry transformation between the current image or frame and the previous image or frame that encodes the rotation and scaled translation information, and R(t) is a representation of the rotation between the current image or frame and the previous image or frame as a function of time. In Equation (3), H can be determined when m_j* and the depth ratio α_j (the ratio of the depth to a point in the current image or frame and the depth to the same point in the previous image or frame) are known Next, the rotation matrix R(t) is decomposed into the angle-axis relationship defined by Equation (12), as indicated in block 124. This transforms R(t) from a 3×3 matrix representation to an angle-axis representation. At this point, e_ω is known (from Equation (12)) and the {circumflex over (ė)}_ω vector can be computed using Equation (14) by selecting gain matrices K_ω and ρ_ω, as indicated in block 126. The {circumflex over (ė)}_ω, vector, which is an estimate of the derivative of e_ω, is determined to filter out noise created by the camera. The angular velocity ω of the camera can then be computed using Equation (16), as indicated in block 128. In that equation, L_ω is an invertible Jacobian matrix.

[0056] With reference to block 130, an actual estimate for the state vector [ŷ₃, û₁, û₂] is computed using Equation (17) and, with reference to block 132 of FIG. 1C, [ y₃, ū₁, ū₂] are computed by integration of Equation (18). As can be appreciated from the foregoing discussion, Equations (17) and (18) are coupled. The initialized [ y₃, ū₁, ū₂] can be used in Equation (16) along with the second term (which is obtained from the image sensor) in that equation to compute [ŷ₃, û₁, û₂] for the initial frame. Then those values can be input into Equation (18) to obtain the derivative of [ y₃, ū₁, ū₂] for the initial frame. Next, [ y₃, ū₁, ū₂] can then be input back into Equation (17) to obtain [ŷ₃, û₁, û₂] for the current frame. As indicated in block 134, ŷ₃ can then be inverted to obtain the z coordinate of each feature point in the current camera coordinate frame. The relations {circumflex over (b)}₁=û₁/ŷ₃ and {circumflex over (b)}₂=û₂/ŷ₃ can then be used to compute the remaining two unknown camera velocities (e.g., the linear velocities in the x and y directions), as indicated in block 136. At this point, the x and y coordinates of each feature point in the current camera coordinate frame can be computed using X=y₁/ŷ₃Y=y₂/ŷ₃, as indicated in block 138. With that information, the relative positions, and motion, of the feature points at the time the current image or frame was captured are known. As indicated in decision block 140, if the estimation is to continue and a further image is to be captured, flow returns to block 112 of FIG. 1A.

[0057] The above-described methods are useful in estimating the structure and motion of stationary objects, but cannot be used for moving objects. Described below are methods for estimating the 3D structure (with proper scale) of objects, whether they are stationary or moving, viewed by a moving camera. In contrast to traditionally developed batch solutions for this problem, a nonlinear unknown input observer strategy is used where the object's velocity is considered as an unknown input to the perspective dynamical system. The estimator is exponentially stable and hence provides robustness against modeling uncertainties and measurement noise from the camera. In some embodiments, the methods implement a first-causal, observer-based structure estimation algorithm for a moving camera viewing a moving object with unknown time-varying object velocities.

[0058] As noted above, recovering the structure of a stationary object viewed by a moving camera is called structure from motion (SfM). The fundamental concept behind SfM algorithms is triangulation. Since the object being observed is stationary, two rays projected onto consecutive images and the camera baseline form a triangle. If the object is not stationary, however, then the projections of the object will be from different locations in the fixed inertial frame. Therefore, triangulation is not feasible and standard SfM techniques cannot be used to recover the structure of a moving object using a moving camera. The following discussion presents a result for the case when the viewed object is in motion. This problem can be described as a structure and motion from motion (SaMfM) estimation.

[0059] Batch algorithms use an algebraic relationship between 3D coordinates of points in the camera coordinate frame and corresponding two-dimensional (2D) projections on the image frame collected over multiple images to estimate the structure. Therefore, batch algorithms are not useful in real-time control algorithms. For visual servo control or video-based surveillance tasks, online structure estimation algorithms are beneficial. The objective of an unknown input observer (UIO) for SaMfM is to estimate the structure of moving objects from a stream of images described using a continuous dynamical model, instead of algebraic relationships and geometric constraints. As described below, a UIO provides an online algorithm that relaxes the constant object linear velocity assumption and does not require the time-varying object velocities to be approximated by a basis.

[0060] Several UIO algorithms are present in literature for estimating the state when an exogenous time-varying unknown input is present in the system. The contribution here is to provide a causal algorithm for estimating the structure of a moving object using a moving camera with relaxed assumptions on the object's motion. In the relative rigid body motion dynamics, the moving object's linear velocity can be viewed as an exogenous time-varying disturbance and is considered as an unknown input. Given the unique UIO method, no assumptions are made on the minimum number of tracked feature points or minimum number of views required to estimate the structure. In fact, the developed approach only requires single tracked feature point and camera velocity information.

[0061] A UIO algorithm for SaMfM will now be presented. Specifically, objectives are defined, a pinhole camera model is introduced and a differential equation formulation of the problem is provided. A UIO is then designed to estimate the unknown structure given measurements from camera images and camera velocities. A procedure for selecting observer gains is developed using a linear matrix inequality (LMI) formulation. Finally, numerical simulations demonstrate the effectiveness of the algorithm.

[0062] The objective of SaMfM is to recover the structure (i.e., Euclidean coordinates with a scaling factor) and motion (i.e., velocities) of moving objects observed by a moving camera, when all camera velocities are assumed to be known. An observer is presented below that estimates the structure of a moving object with respect to an intrinsically calibrated moving camera. At least one feature point on the object is assumed to be tracked in each image frame, and camera velocities are recorded using sensors such as an Inertial Measurement Unit (IMU).

[0063] Consider the scenario depicted in FIG. 2 where a moving camera views moving point objects (such as feature points on a rigid object). In FIG. 2, an inertial reference frame is denoted by F*. The reference frame F* can be attached to the camera at the location corresponding to an initial point in time t where the object is in the camera field of view (FOV), i.e., F* is identical to F_c(t₀). After the initial time, a reference frame F_c attached to a pinhole camera undergoes some rotation R(t)εSO(3) and translation x_ƒ(t)ε³ away from F*.

[0064] The Euclidean coordinates m(t)ε³ of a point observed by a camera expressed in the camera frame F_c and the respective normalized Euclidean coordinates m(t)ε³ are defined as

$\begin{array}{cc}\stackrel{\_}{m}\left(t\right)\stackrel{\Delta }{=}{\left[\begin{array}{ccc}X\left(t\right),& Y\left(t\right),& Z\left(t\right)\end{array}\right]}^T,& \left(27\right)\\ m\left(t\right)\stackrel{\Delta }{=}{\left[\begin{array}{ccc}\frac{X\left(t\right)}{Z\left(t\right)},& \frac{Y\left(t\right)}{Z\left(t\right)},& 1\end{array}\right]}^T.& \left(28\right)\end{array}$

[0065] Consider a closed and bounded set Y⊂³. To facilitate the subsequent development, the state vector x(t)=[x₁(t), x₂ (t), x₃ (t)]^TεY is constructed from Equation (28) as

$\begin{array}{cc}x={\left[\begin{array}{ccc}\frac{X}{Z},& \frac{Y}{Z},& \frac{1}{Z}\end{array}\right]}^T.& \left(29\right)\end{array}$

Using projective geometry, the normalized Euclidean coordinates m(t) can be related to the pixel coordinates in the image space as

q=A_cm  (30)

where q(t)=[ū(t) ν(t) 1]^T is a vector of the image-space feature point coordinates ū(t), ν(t)ε defined on the closed and bounded set I⊂³, and A_cε^3×3 is a constant, known, invertible intrinsic camera calibration matrix. Since A_c is known, the expression in Equation (30) can be used to recover m(t), which can be used to partially reconstruct the state x(t) so that the first two components of x(t) can be determined. The states defined in Equation (29) contain unknown structure information about the object. The object structure can be determined by estimating the Euclidean coordinates of multiple feature points on the object; however, without loss of generality the following UIO development is applied to a single feature point, with Euclidean coordinates that can be recovered from the state x(t) in Equation (29).

[0066] Assumption 5:

[0067] The relative Euclidean distance Z(t) between the camera and the feature points observed on the object is upper and lower bounded by some known positive constants (i.e., the object remains some finite distance away from the camera). Therefore, the definition in Equation (29) can be used to assume that

x₃≧x₃(t)≧x₃  (31)

where x₃, x₃ε denote known positive bounding constants. Likewise, since the image coordinates are constrained (i.e., the object is assumed to remain within the camera field of view (FOV)), the relationships in Equations (28)-(30) along with the fact that A_c is invertible can be used to conclude that

x₁≧|x₁(t)|≧x₁, x₂≧|x₂(t)|≧x₂

where x₁, x₂, x₁, x₂ε denote known positive bounding constants.

[0068] Consider a moving camera viewing a moving point q. As shown in FIG. 2, the point q can be expressed in the coordinate system F_c as

m= x_ƒ+ Rx_oq  (32)

where x_oq is a vector from the origin of coordinate system F* to the point q expressed in the coordinate system F*. Differentiating Equation (32), the relative motion of q as observed in the camera coordinate system can be expressed by the following kinematics

{dot over (m)}=[ω]_x m+ν_r  (33)

where m(t) is defined in Equation (27), [ω]_xε^3×3 denotes a skew symmetric matrix formed from the angular velocity vector of the camera ω(t)=[ω₁ ω₂ ω₃]^TεW, and ν_r(t) represents the relative velocity of the camera with respect to the moving point, defined as

$\begin{array}{cc}{v}_r\stackrel{\Delta }{=}v_c-\stackrel{\_}{R}{\stackrel{\_}{v}}_p.& \left(34\right)\end{array}$

[0069] In Equation (34), ν_c(t)=[ν_cx ν_cy ν_cz]^TεV_c⊂³ denotes the camera velocity expressed in camera reference frame F_c, and

$\stackrel{\_}{R}{\stackrel{\_}{v}}_p\left(t\right)\stackrel{\Delta }{=}v_p\left(t\right)={\left[\begin{array}{ccc}{v}_{\mathrm{px}}& v_{\mathrm{py}}& v_{\mathrm{pz}}\end{array}\right]}^T\in V_p⋐{ℝ}^3$

denotes the velocity of the moving point q expressed in camera reference frame F_c, and

${\stackrel{\_}{v}}_p\left(t\right)\stackrel{\Delta }{=}{\left[\begin{array}{ccc}{\stackrel{\_}{v}}_{\mathrm{px}}& {\stackrel{\_}{v}}_{\mathrm{py}}& {\stackrel{\_}{v}}_{\mathrm{pz}}\end{array}\right]}^T\in {\stackrel{\_}{V}}_p⋐{ℝ}^3$

denotes the velocity of the moving point in the inertial reference frame F*. The sets W, V_c, and V_p are closed and bounded sets such that W⊂³, V_c⊂³, and V_p⊂³. Let the linear and angular camera velocities be denoted by u=[ν_c ω]^T.

[0070] The states defined in Equation (29) contain unknown structure information of the object. To recover the 3D structure, the state x(t) should be estimated. Using Equations (29) and (33), the dynamics of the state vector x(t) are expressed as

{dot over (x)}₁=Ω₁+(ν_cx−x₁ν_cz)x₃−ν_pxx₃+x₁ν_pzx₃,

{dot over (x)}₂=Ω₂+(ν_cy−x₂ν_cz)x₃−ν_pyx₃+x₂ν_pzx₃,

{dot over (x)}₃=−ν_czx₃²−(x₂ω₁−x₁ω₂)x₃+ν_pzx₃²,

y=[x₁x₂]^T.  (35)

where Ω₁(u,y), Ω₂ (u,y)ε are defined as

${\Omega }_1\left(u,y\right)\stackrel{\Delta }{=}-x_1x_2{\omega }_1+{\omega }_2+x_1^2{\omega }_2-x_2{\omega }_3,{\Omega }_2\left(u,y\right)\stackrel{\Delta }{=}-{\omega }_1-x_2^2{\omega }_1+x_1x_2{\omega }_2+x_1{\omega }_3,$

[0071] Assumption 6:

[0072] The camera velocities ω(t), and ν_c(t), the object velocity ν_p(t), and the feature points y(t) are assumed to be upper bounded by constants.

[0073] Remark 5:

[0074] The dynamics in Equation (35) are not observable if: a) the camera is stationary, i.e., ν_c(t)=0,ω(t)=0, or b) the camera moves along the ray projected by the feature point on the image, i.e., (ν_cx−ν_px−x₁(ν_cz−ν_pz))=(ν_cy−ν_py−x₂(ν_cz−ν_pz))=0.

[0075] In the discussion that follows, a UIO is developed for a class of nonlinear systems. Additionally, conditions on the moving object's velocity are stated so that the model in Equation (35) satisfies the requirements of the UIO. The system in Equation (35) can be represented in the form

{dot over (x)}=ƒ(x,u)+g(y,u)+Dd

y=Cx  (36)

where x(t)ε³ is a state of the system, u(t)ε⁶ is a measurable control input, d(t)ε is an unmeasurable disturbance input, y(t)ε² is output of the system, the functions ƒ:³×⁶→³ and g:²×⁶→³ are given by

$f\left(x,u\right)=\left[\begin{array}{c}\left(v_{\mathrm{cx}}-x_1v_{\mathrm{cz}}\right)x_3\\ \left(v_{\mathrm{cy}}-x_2v_{\mathrm{cz}}\right)x_3\\ -x_3^2v_{\mathrm{cz}}-x_2x_3{\omega }_1+x_1x_3{\omega }_2\end{array}\right],g\left(y,u\right)=\left[\begin{array}{c}-x_1x_2{\omega }_1+\left(1+x_1^2\right){\omega }_2-x_2{\omega }_3\\ -\left(1+x_2^2\right){\omega }_1+x_1x_2{\omega }_2+x_1{\omega }_3\\ 0\end{array}\right],C=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right]$

is full row rank and Dε^3×1 is full column rank. The function ƒ(x,u) satisfies following Lipschitz conditions

∥ƒ(x,u)−ƒ({circumflex over (x)},u)∥≦γ₁∥x−{circumflex over (x)}∥  (37)

∥ƒ(x,u)−ƒ({circumflex over (x)},u)−A(x−{circumflex over (x)})∥≦(γ₁+γ₂)∥x−{circumflex over (x)}∥  (38)

where γ₁, γ₂ε⁺, Aε^3×3, and {circumflex over (x)}(t)ε³ is an estimate of the unknown state x(t). Larger γ₁ would mean that the UIO is stable even in the presence of fast moving nonlinear dynamics in ƒ(x,u), i.e., for larger camera velocities.

[0076] The UIO objective is to design an exponentially converging state observer to estimate x(t) in the presence of an unknown input d(t). To quantify this objective an estimation error is defined as

$\begin{array}{cc}\left(t\right)\stackrel{\Delta }{=}\hat{x}\left(t\right)-x\left(t\right).& \left(39\right)\end{array}$

Based on Equation (39) and the subsequent stability analysis, the UIO for the system in Equation (36) is designed as

ż=Nz+Ly+M ƒ({circumflex over (x)},u)+Mg(y,u)

{circumflex over (x)}=z−Ey  (40)

where ƒ(x,u)=ƒ(x,u)−Ax, z(t)ε³ is an auxiliary signal, the matrices Nε^3×3, Lε^3×2, Eε^3×2, Mε^3×3 are designed as

M=I₃+EC

N=MA−KC

L=K(I₂+CE)−MAE  (41)

where E is subsequently designed and Kε^3×2 is a gain matrix which satisfies the inequality

Q=N^TP+PN+(γ₁+γ₂)²PMM^TP+2I₃2I₃<0_3×3  (42)

where Pε^3×3 is a positive definite, symmetric matrix. Using Equation (41), the equality

NM+LC−MA=0_3×3  (43)

is satisfied, where 0_3×3 denotes a zero matrix of the dimensions 3×3. E is selected as

E=F+YG  (44)

where Yε^3×2 can be chosen arbitrarily, and F and G are given by

F=−D(CD)^†,G=(I₂−(CD)(CD)^†).

From the definitions of E, F, and G, ECD=−D and the following equality is satisfied:

MD=(I₃+EC)D=0_3×1.  (45)

Since rank(CD)=1, the generalized pseudo inverse of the matrix CD exists and is given by

(CD)^†=((CD)^T(CD))⁻¹(CD)^T.

Substituting Equation (40) into Equation (39), taking the time derivative of the result, and using Equations (36) and (40) yields

ė=Nz+Ly+M ƒ({circumflex over (x)},u)−(I₃+EC)Ax−(I₃+EC) ƒ(x,u)−(I₃+EC)Dd.  (46)

Using Equations (39) and (41), and the conditions in Equations (43) and (45) the error system in Equation (46) can be written as

ė=Ne+M( ƒ({circumflex over (x)},u)− ƒ(x,u))  (47)

[0077] Theorem 2:

[0078] The nonlinear unknown input observer in Equation (40) is exponentially stable if the assumptions 1 and 2, and condition in Equation (42) is satisfied.

[0079] Proof:

[0080] Consider a positive definite (PD), radially unbounded Lyapunov candidate function V(e): ⁿ→ defined as

V=e^TPe.  (48)

which satisfies

λ_min(P)∥e∥²≦V≦λ_max(P)∥e∥²  (49)

where λ_min(∩) and λ_max(•) are the minimum and maximum eigenvalues of the matrix P. Taking the time derivative of Equation (48) along the trajectories of Equation (47) yields

{dot over (V)}=e^T(N^TP+PN)e+2e^TPM( ƒ({circumflex over (x)},u)− ƒ(x,u))

{dot over (V)}≦e^T(N^TP+PN)e+2∥e^TPM∥∥ƒ({circumflex over (x)},u)−ƒ(x,u)−A({circumflex over (x)}−x)∥

{dot over (V)}≦e^T(N^TP+PN)e+2∥e^TPM∥γ₁∥e∥+2∥e^TPM∥γ₂∥e∥.

Using the norm inequality

2γ_i∥e^TPM∥∥e∥≦γ_i²∥e^TPM∥²+e∥²,∀i={1,2},

the upper bound on {dot over (V)} is given by

{dot over (V)}≦e^T(N^TP+PN)e+(γ₁²+γ₂²)e^TPMM^TPe+2e^Te

{dot over (V)}≦e^T(N^TP+PN+(γ₁²+γ₂²)PMM^TP+2I₃)e

{dot over (V)}≦e^TQe.  (50)

[0081] If the condition in Equation (42) is satisfied, then {dot over (V)}<0. Using Equations (48)-(50) the upper bounds for V(e) and V(e) can be developed as

V(e)≦(e(0))exp(−ξt)

where

$\xi =\frac{{\lambda }_{\max }\left(Q\right)}{{\lambda }_{\min }\left(P\right)}$

and

∥e(t)∥≦ζ∥e(0)∥exp(−ξt)

where

$\zeta =\frac{{\lambda }_{\max }\left(P\right)}{{\lambda }_{\min }\left(P\right)}.$

[0082] Remark 6:

[0083] Model uncertainties can be represented by an additive disturbance term d₁(t)ε³ in Equation (36). Also, if the number of outputs, p, is less than the number of unknown disturbance inputs, q, then q₁<p number of unknown inputs can be represented by d(t) and remaining q−q₁ can be represented by an additive disturbance D₁d₁(t)ε^q−q1. In each of the above mentioned cases, the estimation error will be uniformly ultimately bounded.

[0084] Conditions will next be developed for which an exact estimation of the 3D position of the feature point can be found. Conditions for choosing matrix A in Equation (41) are stated and conditions on the object's trajectories are developed so that an exact estimation can be obtained. The UIO can be used even if the conditions on the object's velocity are not satisfied. In this special case, the state x(t) can only be estimated up to a bounded estimation error, which can be reduced using an optimization procedure.

[0085] Next, conditions for selecting A are provided based on the observer existence conditions. Sufficient conditions for the existence of the UIO in Equation (40) can be summarized as: (MA, C) is observable, rank(CD)=rank(D)=q, MA-KC is Hurwitz, and Equation (42) holds. Since rank(CD)=rank(D)=1, subsequent development indicates that observability of the pair (MA,C) is equivalent to the following rank condition

$\begin{array}{cc}\mathrm{rank}\left[\begin{array}{cc}\lambda I_3-A& D\\ C& 0_{2\times 1}\end{array}\right]=4,\forall \lambda \in C.& \left(51\right)\end{array}$

[0086] Thus, A in Equation (41) should be selected so that Equation (51) is satisfied and (MA,C) is observable. The condition in Equation (51) facilitates the selection of A based on the system matrices and hence circumvents the computation of M for checking the observability of(MA,C). Another criteria on the selection of A is to minimize the Lipschitz constant in Equation (37).

[0087] Some assumptions on the motion of the camera and the viewed object are required to satisfy the constraints on disturbance inputs in Equation (36). In this section, two specific scenarios are discussed.

[0088] The camera is undergoing arbitrary purely translational motion, i.e., the angular velocities of the camera are zero and the object is moving along a straight line with time-varying unknown velocities. For this case, choices of R(t) and ν_p(t) in Equation (34) become R=I₃, and ν_p(t)=[d₁(t) 0 0]^T, or ν_p (t)=[0 d₁(t) 0]^T, or ν_p(t)=[0 0 d₁(t)]^T, or ν_p(t)=[d₁(t) d₁(t) 0]^T, etc., where d₁(t)ε is the unknown time-varying object velocity. Physically, these object velocities would mean the object is moving along straight lines or moving along a circle.

[0089] A downward looking camera is undergoing arbitrary translational motion along with the angular velocity along the Z-direction (an axis pointing downwards). The object is moving on a ground plane (i.e., X-Y plane) with arbitrary time-varying velocities. In this case, the rotation matrix R(t) is given by

$\stackrel{\_}{R}=\left[\begin{array}{ccc}\cos \left(\theta \left(t\right)\right)& \sin \left(\theta \left(t\right)\right)& 0\\ -\sin \left(\theta \left(t\right)\right)& \cos \left(\theta \left(t\right)\right)& 0\\ 0& 0& 1\end{array}\right]=\left[\begin{array}{cc}{\stackrel{\_}{R}}_s& 0_{2\times 1}\\ 0_{1\times 2}& 1\end{array}\right]$

where θ(t)ε[−2π, 2π) is the rotation angle between the current camera coordinate frame and inertial coordinate frame, and R_s(t)ε^2×2 represents the upper left 2×2 block of the R(t). The object velocity in the inertial frame is represented as ν_p(t)=[ d₁(t) d₂(t) 0]^T, where d₁(t), d₂(t)ε. The camera angular velocity is such that

Rν_p=ν_p  (52)

where ν_p1(t)=[ d₃(t) 0 0]^T or ν_p2(t)=[0 d₄(t) 0]^T, and d₃(t), d₄(t)ε are unknown time-varying quantities. The equality in Equation (52) can be achieved if R_s[ d₁(t) d₂(t)]^T=ν_p3, where ν_p3(t)=[ d₃(t) 0]^T or ν_p3(t)=[0 d₄(t)]^T. For example, consider an object undergoing a circular motion with unknown and possibly time-varying radius in the X-Y plane with unknown time-varying velocities observed by a camera undergoing an arbitrary linear motion in the X-Y-Z plane and circular motion along the downward looking Z-direction.

[0090] Next, the condition in Equation (42) is reformulated as an LMI feasibility problem. The matrices P, K and Y should be selected such that the sufficient condition for the observer error stability in Equation (42) is satisfied. Substituting N and M from Equation (41) into Equation (42) yields

(MA−KC)^TP+P(MA−KC)+2I₃+(γ₁²+γ₂²)P(I₃+EC)(I₃+EC)^TP<0_3×3.  (53)

After using Equation (44), the inequality in Equation (53) can be expressed as

A^T(I₃+FC)^TP+P(I₃+FC)A+A^TC^TG^TP_Y^T+P_YGCA−C^TP_K^T−P_KC+2I₃+(γ₁²+γ₂²)(P+PFC+P_YGC)(P+PFC+P_YGC)^T<0_3×3  (54)

where P_Y=PY and P_K=PK. For the observer synthesis, the matrices Y, K and P>0 should be computed such that the matrix inequality in Equation (42) is satisfied. Since P>0, exists, and Y and K can be computed using Y=P⁻¹P_Y, and K=P⁻¹P_K. Using Schur's complement, the inequality in Equation (54) can be transformed into the matrix inequality

$\begin{array}{cc}\left[\begin{array}{cc}{P}_1& \beta R\\ \beta R^T& -I_3\end{array}\right]<0_{6\times 6}& \left(55\right)\end{array}$

where

P₁=A^T(I₃+FC)^TP+P(I₃+FC)A+A^TC^TG^TP_Y^T+P_YGCA−C^TP_K^T−P_KC+2I₃,R=P+PFC+P_YGC,

β=√{square root over (γ)}₁²+γ₂².

[0091] The matrix inequality in Equation (55) is an LMI in variables P, P_Y, and P_K. The LMI feasibility problem can be solved using standard LMI algorithms and is a problem of finding P, P_Y, and P_K such that β is maximized. Maximizing β is equivalent to maximizing γ, which means the observer can be designed for nonlinear functions with a larger Lipschitz constant. If the LMI in Equation (55) is feasible, then a solution to Equation (42) exists. Hence, the LMI feasibility problem is a sufficient condition for the stability of the observer.

[0092] Remark 7:

[0093] In Equation (40), the gain K is multiplied by the measurement y(t). To reduce the effects of noise, the elements of gain K should not be large. Since K=P⁻¹P_K, the minimum eigenvalue of P should be maximized and the maximum eigenvalue of P_K should be minimized. This can be achieved via a multi-objective optimization problem with the objective function as maximize δλ_min(P)−(1−δ)λ_max(P_K) subject to the LMI in Equation (55), where λε[0,1]. In the presence of model uncertainties or disturbance inputs as stated in Remark 2, another optimization objective can be to choose gain K such that the LL₂ gain from the disturbance input to the estimation error is minimized.

[0094] Two simulations are performed for a moving camera observing an object moving in a plane. For the first simulation, camera velocities are given by ν_c(t)=[2 1 0.5 cos(t/2)]^T m/s, and ω(t)=[0 0 1]^T rad/s. The object velocity is selected such that ν_p(t)=[0.5 sin(t) 0 0]^T m/s. The camera calibration matrix is given by

$A_c=\left[\begin{array}{ccc}720& 0& 320\\ 0& 720& 240\\ 0& 0& 1\end{array}\right].$

Measurement noise with 20 dB SNR is added to the image pixel coordinates and the camera velocities using awgn( ) command in Simulink. Matrices A, C, and D are

$A=\left[\begin{array}{ccc}0& -1& 2\\ 1& 0& 1\\ 0& 0& 0\end{array}\right],C=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right],D=\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right].$

The matrix Y and the gain matrix K are computed using the LMI feasibility command ‘feasp’ in Matlab and are given by

$K=\left[\begin{array}{cc}0.8278& 0\\ 0& 0.8278\\ -1.5374& 0\end{array}\right],Y=\left[\begin{array}{cc}0& 0\\ 0& -1\\ 0& -1.5374\end{array}\right].$

FIG. 3 shows comparison of the actual and estimated X, Y and Z coordinates of the object in the camera coordinate frame.

[0095] In the second simulation, camera velocities are given by ν_c(t)=[2 1+0.02t −0.01]^T m/s, ω(t)=[0 0 0.1]^T rad/s. The object motion is set to a circular motion with unknown radius and velocities. Measurement noise with 20 dB SNR is added to the image pixel coordinates and the camera velocities using awgn( ) command in Simulink. Matrices A, C and D are:

$A=\left[\begin{array}{ccc}0& -0.1& 2\\ 0.1& 0& 1\\ 0& 0& 0\end{array}\right],C=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right],D=\left[\begin{array}{c}0\\ 1\\ 0\end{array}\right].$

The matrix Y and the gain matrix K are computed using the LMI feasibility command ‘ feasp’ in Matlab and are given by

$K=\left[\begin{array}{cc}1.0667& 0\\ 0& 1.0667\\ 0& 0.056\end{array}\right],Y=\left[\begin{array}{cc}-1& 0\\ 0& 0\\ -0.5598& 0\end{array}\right].$

FIG. 4 shows the comparison of actual and estimated X, Y, and Z coordinates of the object in the Inertial coordinate frame in the presence of measurement noise. FIG. 5 shows the structure estimation error. From the two results in two simulations, it can be seen that the observer is insensitive to the disturbance input.

[0096] Experiments are conducted on PUMA 560 serial manipulator to test the performance of the proposed estimator. Results of one of the test cases are presented below. A mvBlueFox-120a USB camera is rigidly attached to the end-effector of PUMA. A two-link robot mounted on the same platform as the PUMA is used as a moving object. The PUMA is controlled using a workstation to generate camera motion, and the camera linear and angular velocities are computed using the joint encoders and forward velocity kinematics of the PUMA. The data from the camera is recorded on the PUMA workstation computer at 25 frames per second for further processing and is time synchronized with the camera images. Two-link robot was controlled using another workstation which was time-synchronized with the PUMA workstation for ground truth verification. The video data recorded from the camera is post-processed to track a feature point on second link of the two-link robot using KLT feature point tracker. The image processing code is implemented in Visual C++ 2008 using OpenCV library. Since the PUMA and two-link robot was mounted on the same table, a precise location of the base of both the robots is computed and forward kinematics is used to compute the position of the camera and the feature point on a link of the two-link robot. Hence, a 3D distance between camera mounted on the PUMA and feature point on the two-link robot can be computed precisely and is used as a ground truth for testing the performance of the proposed estimator.

[0097] A downward looking camera attached to PUMA is moved with randomly chosen linear and angular velocities in X-Y-Z directions, and the two-link robot is commanded certain velocities in X-Y plane. The observer in Equations (40)-(41) is implemented to estimate the relative 3D location of a moving feature point on the two-link robot in the camera reference frame. The observer parameters are selected as

$\begin{array}{cc}A=\left[\begin{array}{ccc}0& -0.05& 0\\ 0.05& 0& -0.3\\ 0& 0& 0\end{array}\right],C=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right],D=\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right].& \left(56\right)\end{array}$

The matrix Y and the gain matrix K are computed using the LMI feasibility command feasp in Matlab to solve LMI in Equation (55) using the matrices in Equation (56) and are given by

$K=\left[\begin{array}{cc}1.2204& 0\\ 0& 1.2204\\ 0.3731& 0.056\end{array}\right],Y=\left[\begin{array}{cc}0& 0\\ 0& -1\\ 0& 7.4618\end{array}\right].$

FIGS. 6 and 7 shows the camera linear and angular velocities. The comparison of the estimated relative 3D position of the moving point in the camera reference frame and ground truth 3D position is shown in FIG. 8.

[0098] Some preliminary studies of the structure (3D location) of a moving object using a moving camera with known camera velocities has been presented. More testing and experimental verification of the proposed algorithm is required to verify the performance of the algorithm in various real-life situations.

[0099] FIGS. 9A and 9B illustrate an example method for estimating the structure and motion of a stationary or moving object consistent with the above discussion. In this method, the angular and linear velocities of the camera are known. Beginning with block 200, the auxiliary vector [{circumflex over (z)}₁, {circumflex over (z)}₂, {circumflex over (z)}₃] is initialized. The auxiliary vector is a transformed version of the state vector [{circumflex over (x)}₁, {circumflex over (x)}₂, {circumflex over (x)}₃]. As with the previous method, substantially any initial values (e.g., random numbers) can be used to initialize the vector, although the closer the selected values are to the actual values that are being estimated, the faster convergence will be attained.

[0100] Once the auxiliary vector has been initialized, a matrix D is chosen based on what is believed to be the relative motion between the object and the camera, as indicated in block 202. Matrix D is a representation of that relative motion. Next, with reference to block 204, a matrix A is chosen according to Equation (51). Matrix A is the state matrix and represents camera motion, i.e., the linear and angular velocities of the camera. Notably, the method would work for an arbitrary A matrix also, although it works more efficiently when A represents camera motion. Referring next to blocks 206-210, matrices F and G can be computed using Equation (54), matrices K and Y can be computed using Equation (55), matrices M, N, and L can be computed using Equation (51), and matrix E can be computed using Equation (54), wherein matrices E, F, G, M, N, and L are observer parameters and matrices K and Y are gain matrices.

[0101] Turning to block 212, an image of the object can be captured using a moving camera. Again, FIG. 2 schematically illustrates such image capture with a camera. By way of example, the object is a moving object and the camera is mounted to a further moving object, such as a UAV. Once the image is captured, feature points within the image are extracted, as indicated in block 214. By way of example, four coplanar feature points are extracted or eight non-coplanar feature points are extracted. In some embodiments, coplanar points can be identified using image segmentation techniques. With reference to block 216, the x and y pixel coordinates of each image feature point are recorded and, with reference to block 218, the linear velocity and linear acceleration of the camera motion at the time of image capture are recorded.

[0102] At this point, the Euclidean coordinates m_j of the feature points are computed, as indicated in block 220. By way of example, those coordinates are determined using Equation (30). Referring next to block 222, an estimate for the state vector [{circumflex over (x)}₁, {circumflex over (x)}₂, {circumflex over (x)}₃] is computed using the second relation of Equation (40), which relates the state vector to the pixel coordinates y, the auxiliary vector z, and the matrix E. As can be appreciated from the foregoing discussion, {circumflex over (x)}₁, is the inverse depth of the feature point (i.e., 1/Z), {circumflex over (x)}₂ is a ratio of the x and z coordinates (i.e., X/Z), and {circumflex over (x)}₃, is a ratio of the y and z coordinates (i.e., Y/Z. Next, X, is inverted to obtain the z coordinate of each feature point in the current camera coordinate frame, as indicated in block 224, and the relations X=x₁/{circumflex over (x)}₃ and Y=x₂/{circumflex over (x)}₃, are used to compute the x and y coordinates of each feature point in the current camera coordinate frame, as indicated in block 226.

[0103] At this point, auxiliary vector [{circumflex over (z)}₁,{circumflex over (z)}₂,{circumflex over (z)}₃] is integrated using matrices N, L, M, and E, as indicated in block 228, to obtain a new auxiliary vector that can be used in the next iteration of the process (beginning with block 212). As indicated in decision block 230, if the estimation is to continue and a further image is to be captured, flow returns to block 212 of FIG. 9A.

[0104] FIG. 10 illustrates an example architecture for a computing device 300 that can be used to perform at least part of the processes described above. As indicated in FIG. 10, the computing device 300 at least comprises a processing device 302 and memory 304. The processing device 302 can include a central processing unit (CPU) or other processing device (e.g., microprocessor or digital signal processor) and the memory 304 includes any one of or a combination of volatile memory elements (e.g., RAM) and nonvolatile memory elements (e.g., flash, hard disk, ROM).

[0105] The memory 304 is a non-transitory computer-readable medium that stores various programs (i.e., logic), including an operating system 306 and an structure and motion estimator 308. The operating system 306 controls the execution of other programs and provides scheduling, input-output control, file and data management, memory management, and communication control and related services. The structure and motion estimator 308 comprises one or more algorithms that are configured to receive image data and analyze that image data to estimate the structure and motion of objects captured in the image data.

[0106] Although particular embodiments of systems and methods have been described in this disclosure, it is to be understood that various other embodiments are possible. All such embodiments are intended to fall within the scope of this disclosure.