Ellipse Fitting Based Approach for Extended Object Tracking

With the increase of sensors’ resolution, traditional object tracking technology, which ignores object’s physical extension, gradually becomes inappropriate. Extended object tracking (EOT) technology is able to obtain more information about the object through jointly estimating both centroid’s dynamic state and physical extension of the object. Randommatrix based approach is a promising method for EOT. It uses ellipse/ellipsoid to describe the physical extension of the object. In order to reduce the physical extension estimation error when object maneuvers, the relationship between ellipse/ellipsoid and symmetrical positive definite matrix is analyzed at first. On this basis, ellipse/ellipsoid fitting based approach (EFA) for EOT is proposed based on themeasurement model and centroid’s dynamic model of random matrix based EOT approach. Simulation results show that EFA is effective. The physical extension estimation error of EFA is lower than those of random matrix based approaches when object maneuvers. Besides, the estimation error of centroid’s dynamic state of EFA is also lower.


Introduction
In recent years, object tracking has always been the research focus in information fusion and many other fields [1,2].Physical extension of targets can be ignored if sensor's resolution is relatively low.Traditional object tracking techniques usually ignore target's physical extension and use centroid's dynamic state to depict that of the entire target [3].As the resolution of measurement sensors is gradually improved, this hypothesis becomes no longer suitable [4].For example, when tracking large surface ships (such as aircraft carriers), each scan may produce several measurement data and the number of these data varies all the time.When some objects compose an intensive formation, such as flight formation in aerial refueling, the measurement sensors may fail to distinguish every individual object.These situations make it quite difficult to establish one-to-one correspondence between an object and a measurement.Traditional OT technology is no longer appropriate in these cases.Thus, it is necessary to measure and estimate object's physical extension by extended object tracking (EOT) technique to improve the tracking performance [4][5][6][7][8].Different from traditional OT, EOT estimates object's physical extension and centroid's dynamic state jointly.According to different modeling methods of the physical extension, EOT approaches include particle filter [9], probability hypothesis density filter [10,11], and the approach based on random hypersurface model [12].Compared with the above methods, random matrix based EOT approach (RMF) is simpler and has a smaller amount of computation, which makes it more practical and promising.RMF is first proposed by Koch [4] and uses ellipse/ellipsoid to describe the physical extension of object.This simplification can be applied in many scenarios, especially in the military fields, as shown in Figure 1.
Because ellipse/ellipsoid can be expressed by symmetrical positive definite (SPD) matrix, the extension of object can be further depicted by SPD matrix.RMF assumes that the physical extension evolution is subject to Wishart-related distributions [4,7].The basic principle of RMF is to convert the joint posteriori probability distribution to the product of a Gaussian distribution and an inverse Wishart distribution within Bayesian framework.The filtering equations are obtained during the derivation.Feldmann et al. [7] improved measurement noise model and applied RMF in interactive multiple-model algorithm.Lan and Li [5,6] further improved both of physical extension evolution model and measurement noise model.Li et al. [13] proposed the adaptive RMF based on model parameter adaption.According to the simulation results in literature, physical extension estimation error of RMF significantly increased when object is maneuvering, even utilizing multimodel algorithm [4][5][6][7]13].The "maneuvering" here refers to object's actions that influence the physical extension, such as turning [4][5][6][7]13] and significant change of object size.
In order to improve EOT performance and reduce the physical extension estimation error during object's maneuvering, the relationship between ellipse/ellipsoid and SPD matrix is analyzed in this paper.Then ellipse/ellipsoid fitting based approach (EFA) for EOT is proposed.EFA can be integrated into a switching scheme, as shown in [8].Utilizing "Givens rotation" [14] for SPD matrix decomposition, EFA divides the ellipse/ellipsoid fitting problem into two subproblems: direction angle estimation and semiaxis lengths estimation.The simulation results show that either the object is maneuvering or not; physical extension estimation error of EFA does not change much.When the object maneuvers, physical extension estimation error of EFA is lower than RMF.
The object state is characterized by centroid's dynamic state x  and physical extension X  simultaneously in EOT.Their estimation processes in EFA are mutually influenced and will be introduced successively in the sequel.

Centroid's Dynamic State Estimation
Let  denote the spatial dimension, which is usually set as 2 or 3.  is the number of dimensions of centroid's dynamic state in a spatial dimension.For example, when  = 2, the centroid's dynamic state includes position and velocity; when  = 3, it also includes the acceleration [4][5][6][7]13].Thus, the centroid's dynamic state x  is a random  × 1 vector, while the physical extension X  is a random  ×  SPD matrix representing the ellipse/ellipsoid in -dimension space.This section will introduce the dynamic model and measurement model of x  .Then the centroid's dynamic state estimation x | can be obtained by suitable estimators, for example, standard Kalman Filter.

Dynamic Model.
The dynamic model of x  used in EOT is quite similar to that of standard Kalman filter [5][6][7]13] as follows: where F  = F ⊗ I  is the state matrix, F is the state matrix in one spatial dimension, ⊗ is the right Kronecker product [15], and I  is the  ×  identity matrix.  is the independent process noise, which is subject to the following Gaussian distribution: where D is the covariance matrix of   in a spatial dimension and N(, Σ) is the Gaussian distribution with the expectation of  and the covariance matrix of Σ.

Measurement Model.
At scan , the measurement sensor obtains measurement data z   ( = 1, 2, . . .,   ) with   ≥ 1, where   is independent from X  and x  and subject to Poisson distribution with mean  P .The real distribution of these measurement data is set as uniform distribution over the object, which is depicted by ellipse/ellipsoid, with the Gaussian noise of the measurement sensor superimposed [5,7] as follows: Measurement model ( 3) is hard to be directly used in the Bayesian framework or Kalman-related filter.Therefore, the following model of a single measurement is used for approximating [4,7]: where H  = H ⊗I  and H is the measurement matrix in one spatial dimension.Let 0 −1 be ( − 1) × 1 zero vector and H = [1, 0 −1 ]; then, only the position of centroid is measured.
) and is subject to the Gaussian distribution N (0, X  + Ψ  ), where  is a scalar constant. can be set as 0.25 [7] via equaling the covariance matrices of z   in ( 3) and ( 4).Since the probability density of the Gaussian distribution near the expectation is relatively large, only equaling two covariance matrices would lead to excessive measurement data near the ellipse's center.Therefore,  can be set as a relatively larger value so that (4) could be closer to the real measurement model as shown in (3).
Furthermore, the simplified centroid's measurement model for estimation of x  in EFA is as follows: where H  is the same as that in (4) and In (5),   ∼ N(0, X −1|−1 + Ψ  ) is the simplified measurement noise, and X −1|−1 is the estimation of the physical extension at scan  − 1 ( > 1), which is used to approximate X  .Equation ( 5) can be interpreted as taking the mean value of   measurement data as the final measurement result at scan .It indicates that the distribution of measurement data is influenced by object's physical extension X  , which also means that the estimation process of centroid's dynamic state is influenced by physical extension.

Physical Extension Estimation
Since each scan would produce several measurement data in EOT, the basic idea of physical extension estimation in EFA is to treat the position components of these measurement data as points in -dimension space.These measurement points could be used for ellipse/ellipsoid fitting, and then the estimated value of object's physical extension is obtained based on the relationship between ellipse/ellipsoid and SPD matrix.
In this section, EFA firstly uses "Givens rotations" [14] to decompose SPD matrix and analyze the relationship between ellipse/ellipsoid and SPD matrix.On this basis, ellipse/ellipsoid fitting problem is divided into the estimation of direction angle and semiaxis lengths, relatively.
As most cases in EOT research set  as 2 [4-7, 13], we take  = 2 as the example for introduction.The principle of EFA when  = 3 is analogous to that when  = 2.

Ellipse/Ellipsoid and SPD Matrix.
Based on the dynamic model and measurement model of x  defined in Section 2, the equation of ellipse/ellipsoid which depicts the object is as follows [5]: where  × 1 vector ⃗ x is the coordinates of points on the ellipse/ellipsoid.
To facilitate the analysis of the relationship between ellipse/ellipsoid and SPD matrix, the following lemma is introduced at first.Lemma 1. × SPD matrix A can be converted into a diagonal matrix by "Givens rotations" with the maximal times of ( − 1)/2.The diagonal elements of the diagonal matrix are the eigenvalues of A [14].
When  = 2, according to Lemma 1, X  could be decomposed into where are eigenvalues of X  , and An SPD matrix can be decomposed into the form shown in ( 8) by singular value decomposition, eigenvalue decomposition, and other means.However, these means could only ensure that R(  ) is an orthogonal matrix.By Givens rotations, it could be guaranteed that R(  ) is a rotation matrix.If the determinant of an orthogonal matrix is +1, then it is a rotation matrix.
With (7) in mind, the following three elements are needed to completely identify an ellipse/ellipsoid in -dimension space: (a) coordinates of the center point, which are identified by position component of x  , (b) direction angle, which is identified by the rotation angle of "Givens rotations" in Lemma 1: ( − 1)/2 direction angles are needed at most to describe an ellipse/ellipsoid in -dimension space, (c) the semiaxis lengths are equal to the square roots of X  's eigenvalues.
Therefore, X  and x  could jointly depict any ellipse/ellipsoid in -dimension space.When the direction angles and semiaxis lengths of ellipse/ellipsoid are identified, X  could be uniquely determined according to Lemma 1.

Direction Angle Estimation.
Considering that X  could depict all of the ellipses/ellipsoids centering at the origin, according to (8), any ellipse X  on a 2-dimension (2D) plane could be obtained by rotation of ellipse depicted by X 0 = diag( 2 ,1 ,  2 ,2 ), where  2 ,1 and  2 ,2 are eigenvalues of X  .X 0 is regarded as the "reference ellipse" in EFA, whose axes are parallel to two coordinate axes, relatively.If we assume that  2  ,1 ≥  2 ,2 > 0, then the direction angle   in (8) is the intersection angle between the major axis of ellipse X  and the -axis, while the counter-clockwise direction is set as the positive direction.Therefore, the method to estimate ellipse's direction angle   when  = 2 is as follows.
(a) Centralize measurement data: (c) Solve the 2 × 1 feature vector ⃗ r z,max corresponding to the larger eigenvalue of Σ * z .Then the intersection angle between ⃗ r z,max and -axis is the estimated value of   as follows: ) .
Step (a) of direction angel estimation indicates that physical extension estimation is influenced by centroid's dynamic state.According to the principle of principal component analysis [16], ⃗ r z,max is the direction of the largest variance of z , *  .In addition, when utilizing the measurement model shown in (3), the direction of X  's major axis could be approximately considered the direction with the greatest changes of z , *  .Therefore, EFA uses ⃗ r z,max as the estimated value of the major axis direction of X  , which is also the estimation of direction angle   .

Semiaxis Lengths Estimation.
The method of estimating semiaxis lengths in EFA is divided into two steps as follows: (a) identifying the upper and lower bounds of semiaxis lengths; (b) selecting a suitable step size that is selected for traversing, and then the semiaxis lengths of an ellipse are estimated according to corresponding criterion.
When  = 2, the semiaxis lengths of ellipse X  are  ,1 ≥  ,2 > 0. Their upper bounds can be determined as follows: (a) z , *  is projected to the directions of two feature vectors of Σ * z ; (b) the maximal value of absolute projected coordinates could be taken as the upper bounds  max ,1 and  max ,2 , respectively.
The lower bounds of semiaxis lengths could be identified by the covariance matrix Ψ  of measurement sensor's noise.Let  2 Ψ, ≥  2 Ψ, > 0 denote the eigenvalues of Ψ  .Then the lower bound of semiaxis lengths could be set as  min ,1 =  Ψ, and  min ,2 =  Ψ, .Once the upper and lower bounds of semiaxis lengths are identified, the method to estimate semiaxis lengths is as follows.
(a) Determine traversal step size according to the required accuracy and operating speed.
(b) Starting from  min ,1 and  min ,2 , each pair of semiaxis lengths candidates  ,1 and  ,2 that meets the conditions of ( 12) is traversed.And the corresponding estimated X | is calculated according to (8) as follows: (c) Calculate  den based on (7), which is the number of measurement points on and inside the ellipse.
(d) If  den /  >  den , the density of measurement points is calculated as  den ∝  den / ,1  ,2 , which is the number of measurement points in the unit area of ellipse, and 0 <  des < 1 is a preset threshold; else,  den is set to an infinitesimal.
(e) If there is no candidate value of semiaxis lengths that meets (12), advance to Step (f); else, return to Step (b).
(f) The candidate value of semiaxis lengths which has the largest corresponding  den is the estimated value σ,1 and σ,2 .
In Step (d),  den could be replaced by the lower bound of the Wilson interval [17] as follows: where  is set as 1.96 under the confidence level of 95% and ξden is the normalization result of all valid  den .When using  den , the improper situation that  des and  ,1  ,2 are all quite small but that  den is large is easy to occur. wilson simultaneously considers  den and  den ; that is, it is required that the density is maximized while the number of measurement points on and inside the ellipse should not be too small.Thus, the aforementioned improper situation could be effectively avoided.When   is large, the effect of  wilson and  den is more or less equivalent.Since   in EOT is usually not very large,  wilson would perform better than  den in EFA, yet resulting in more complex calculations.
To take into account higher precision and faster calculation speed simultaneously, we could use a relatively larger step size to obtain σ,1 and σ,2 , and then repeat the above algorithm with a smaller step size in the interval near σ,1 and σ,2 .Thus, a more efficient estimation result can be achieved through a relatively smaller amount of computation.
After calculating the direction angle and semiaxis lengths, the estimated value X | of X  could be obtained according to (8).

Physical Extension Estimation
When  = 3.When  = 3, according to Lemma 1, we have where R ⋅ ( ,⋅ ) is the rotation matrix having counterclockwise rotation around a certain coordinate axis.For example, by counterclockwise rotation around -axis, it should be Equation (14) shows that maximally three Givens rotations are needed to convert X  to a diagonal matrix when  = 3.Similar to that when  = 2, we may assume that  2 ,1 ≥  2 ,2 ≥  2 ,3 > 0 and adopt the following sequence of rotation matrix: where R  ( ,1 ) is the rotation matrix with the angle  ,1 around -axis at the counterclockwise direction, which is the same with R  ( ,2 ) and R  ( ,3 ).Therefore, the ellipsoid's direction angles  , ( = 1, 2, 3) could be estimated as follows.
(b) Project z , *  to the - 2D plane.The estimation method of direction angle when  = 2 is used to obtain φ,1 .
When  = 3, the semiaxis lengths estimation method with  = 2 is completely suitable.The difference is that  den ∝  den / ,1  ,2  ,3 represents the number of measurement points in the unit volume of the ellipsoid.
Once the direction angles and semiaxis lengths are obtained, the estimated value X | could be calculated according to (15).

Simulation Results and Discussion
4.1.Scenario 1. Scenario 1, as shown in Figure 2, is commonly used in RMF research [5,7,13].The object is an aircraft carrier with the length of about 350 m and the width of about 100 m.An ellipse with the semiaxis lengths of 175 m and 50 m is used to depict its physical extension.The object starts from the origin and performs uniform linear motion along the path.Then it makes three turns as shown in Figure 2. The velocity of object is set as 27 knots (about 50 km/h).The measurement period is set as  = 10s and the covariance matrix of measurement sensor's noise is Ψ  = diag(40 2 , 15 2 ) m 2 .The number of measurement points of each scan is subject to Poisson distribution with the mean value of  p = 25.
The following four EOT approaches are compared through  MC = 600 Monte Carlo simulations: (a) Koch's RMF [4]: single-model RMF approach and time decay constant  is set as 8; (b) Feldmann's RMF-MM [7]: multimodel RMF approach and parameters are set according to the literature [7]; (c) EFA-den:  den is selected for semiaxis lengths estimation, and threshold value  den is set as 0.5; (d) EFA-wilson:  wilson is selected for semiaxis lengths estimation, and threshold value  den is set as 0.5.
To compare the performance of physical extension estimation, root mean square error (RMSE) of physical extension is defined as follows [7]: where X ℎ | is the estimated value of physical extension at scan  in ℎth Monte Carlo simulation.Physical extension RMSE curves in Scenario 1 are shown in Figure 3.When the object is maneuvering (making three turns), physical extension RMSEs of RMFs increase while those of EFAs do not change significantly.Physical extension RMSEs of EFAs are lower than Koch's RMF but higher than Feldmann's RMF-MM when the object is not maneuvering.Physical extension RMSE of EFA-wilson is lower than that of EFA-den.

Mathematical Problems in Engineering
The centroid's position and velocity RMSE curves in Scenario 1 are shown in Figures 4 and 5.The position and velocity RMSEs of two EFAs are nearly the same and lower than two RMFs.

Scenario 2.
In Scenario 2, the trajectory of object is shown in Figure 6.From the origin, the object takes uniform linear motion for some time and then keeps taking uniform circular motion.Four EOT approaches and other simulation parameters for comparison are the same as those in Scenario 1. RMSE curves of physical extension are shown in Figure 7.In the beginning, the object is not maneuvering, and physical extension RMSE of Feldmann's RMF-MM is lower than EFAs.When the object is taking uniform circular motion (maneuvering), the physical extension RMSEs of two EFAs become lower than two RMF approaches.
RMSE curves of centroid's position and velocity are shown in Figures 8 and 9. Similar to Scenario 1, the position and velocity RMSEs of two EFAs are quite close to each other and lower than the RMFs.

Conclusions
Extended object tracking is able to estimate centroid's dynamic state and physical extension of objects jointly and provide more information about the object.Therefore, it has become a research focus of information fusion and maneuvering object tracking in recent years.RMF uses simple models and could obtain iterative filtering equations within the Bayesian framework.Its main drawback is the significant rise of physical extension estimation error during object's maneuvering.In order to overcome the shortcoming of RMF, an EOT approach based on ellipse/ellipsoid fitting is proposed on the basis of the analysis of the internal correlation between ellipse/ellipsoid and SPD matrix.The simulation results show that EFA is insensitive to whether the object is maneuvering or not.Physical extension estimation error of EFA is lower than the two RMFs during object's maneuvering.Centroid's dynamic state estimation error of EFA is lower than RMF.EFA-wilson performs better than EFA-den in physical extension estimation.

Figure 1 :
Figure 1: Describing physical extension of object with an ellipse: aircraft carrier (a) and flight formation (b).