An extended target tracking model with multiple random matrices and unified kinematics
Karl~Granstr\"om% <-this % stops a space
~\IEEEmembership{Member
~IEEE
....lastname \thanks{...
%\boldmath This paper presents a model for tracking of extended targets, where each target is represented by a given number of elliptic subobjects. A gamma Gaussian inverse Wishart implementation is derived, and necessary approximations are suggested to alleviate the data association complexity. A simulation study shows the merits of the model compared to previous work on the topic.
MultipleEllipse.tex
%
%% bare_jrnl.tex
%% V1.3
%% 2007/01/11
%% by Michael Shell
%% see http://www.michaelshell.org/
%% for current contact information.
%%
%% This is a skeleton file demonstrating the use of IEEEtran.cls
%% (requires IEEEtran.cls version 1.7 or later) with an IEEE journal paper.
%%
%% Support sites:
%% http://www.michaelshell.org/tex/ieeetran/
%% http://www.ctan.org/tex-archive/macros/latex/contrib/IEEEtran/
%% and
%% http://www.ieee.org/
%%


% *** Authors should verify (and, if needed, correct) their LaTeX system  ***
% *** with the testflow diagnostic prior to trusting their LaTeX platform ***
% *** with production work. IEEE's font choices can trigger bugs that do  ***
% *** not appear when using other class files.                            ***
% The testflow support page is at:
% http://www.michaelshell.org/tex/testflow/


%%*************************************************************************
%% Legal Notice:
%% This code is offered as-is without any warranty either expressed or
%% implied; without even the implied warranty of MERCHANTABILITY or
%% FITNESS FOR A PARTICULAR PURPOSE!
%% User assumes all risk.
%% In no event shall IEEE or any contributor to this code be liable for
%% any damages or losses, including, but not limited to, incidental,
%% consequential, or any other damages, resulting from the use or misuse
%% of any information contained here.
%%
%% All comments are the opinions of their respective authors and are not
%% necessarily endorsed by the IEEE.
%%
%% This work is distributed under the LaTeX Project Public License (LPPL)
%% ( http://www.latex-project.org/ ) version 1.3, and may be freely used,
%% distributed and modified. A copy of the LPPL, version 1.3, is included
%% in the base LaTeX documentation of all distributions of LaTeX released
%% 2003/12/01 or later.
%% Retain all contribution notices and credits.
%% ** Modified files should be clearly indicated as such, including  **
%% ** renaming them and changing author support contact information. **
%%
%% File list of work: IEEEtran.cls, IEEEtran_HOWTO.pdf, bare_adv.tex,
%%                    bare_conf.tex, bare_jrnl.tex, bare_jrnl_compsoc.tex
%%*************************************************************************

% Note that the a4paper option is mainly intended so that authors in
% countries using A4 can easily print to A4 and see how their papers will
% look in print - the typesetting of the document will not typically be
% affected with changes in paper size (but the bottom and side margins will).
% Use the testflow package mentioned above to verify correct handling of
% both paper sizes by the user's LaTeX system.
%
% Also note that the "draftcls" or "draftclsnofoot", not "draft", option
% should be used if it is desired that the figures are to be displayed in
% draft mode.
%
\documentclass[journal,twoside]{IEEEtran}
%\documentclass[9pt,technote]{IEEEtran}
%\documentclass[draftcls, onecolumn, 11pt]{IEEEtran}
%
% If IEEEtran.cls has not been installed into the LaTeX system files,
% manually specify the path to it like:
% \documentclass[journal]{../sty/IEEEtran}





% Some very useful LaTeX packages include:
% (uncomment the ones you want to load)


% *** MISC UTILITY PACKAGES ***
%
%\usepackage{ifpdf}
% Heiko Oberdiek's ifpdf.sty is very useful if you need conditional
% compilation based on whether the output is pdf or dvi.
% usage:
% \ifpdf
%   % pdf code
% \else
%   % dvi code
% \fi
% The latest version of ifpdf.sty can be obtained from:
% http://www.ctan.org/tex-archive/macros/latex/contrib/oberdiek/
% Also, note that IEEEtran.cls V1.7 and later provides a builtin
% \ifCLASSINFOpdf conditional that works the same way.
% When switching from latex to pdflatex and vice-versa, the compiler may
% have to be run twice to clear warning/error messages.






% *** CITATION PACKAGES ***
%
\usepackage{cite}
% cite.sty was written by Donald Arseneau
% V1.6 and later of IEEEtran pre-defines the format of the cite.sty package
% \cite{} output to follow that of IEEE. Loading the cite package will
% result in citation numbers being automatically sorted and properly
% "compressed/ranged". e.g., [1], [9], [2], [7], [5], [6] without using
% cite.sty will become [1], [2], [5]--[7], [9] using cite.sty. cite.sty's
% \cite will automatically add leading space, if needed. Use cite.sty's
% noadjust option (cite.sty V3.8 and later) if you want to turn this off.
% cite.sty is already installed on most LaTeX systems. Be sure and use
% version 4.0 (2003-05-27) and later if using hyperref.sty. cite.sty does
% not currently provide for hyperlinked citations.
% The latest version can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/cite/
% The documentation is contained in the cite.sty file itself.






% *** GRAPHICS RELATED PACKAGES ***
%
\ifCLASSINFOpdf
  \usepackage[pdftex]{graphicx}
  % declare the path(s) where your graphic files are
   \graphicspath{}
  % and their extensions so you won't have to specify these with
  % every instance of \includegraphics
  % \DeclareGraphicsExtensions{.pdf,.jpeg,.png}
\else
  % or other class option (dvipsone, dvipdf, if not using dvips). graphicx
  % will default to the driver specified in the system graphics.cfg if no
  % driver is specified.
   \usepackage[dvips]{graphicx}
  % declare the path(s) where your graphic files are
   \graphicspath{}
  % and their extensions so you won't have to specify these with
  % every instance of \includegraphics
   \DeclareGraphicsExtensions{.eps}
\fi
% graphicx was written by David Carlisle and Sebastian Rahtz. It is
% required if you want graphics, photos, etc. graphicx.sty is already
% installed on most LaTeX systems. The latest version and documentation can
% be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/required/graphics/
% Another good source of documentation is "Using Imported Graphics in
% LaTeX2e" by Keith Reckdahl which can be found as epslatex.ps or
% epslatex.pdf at: http://www.ctan.org/tex-archive/info/
%
% latex, and pdflatex in dvi mode, support graphics in encapsulated
% postscript (.eps) format. pdflatex in pdf mode supports graphics
% in .pdf, .jpeg, .png and .mps (metapost) formats. Users should ensure
% that all non-photo figures use a vector format (.eps, .pdf, .mps) and
% not a bitmapped formats (.jpeg, .png). IEEE frowns on bitmapped formats
% which can result in "jaggedy"/blurry rendering of lines and letters as
% well as large increases in file sizes.
%
% You can find documentation about the pdfTeX application at:
% http://www.tug.org/applications/pdftex





% *** MATH PACKAGES ***
%
\usepackage[cmex10]{amsmath}
\usepackage{amssymb}

% A popular package from the American Mathematical Society that provides
% many useful and powerful commands for dealing with mathematics. If using
% it, be sure to load this package with the cmex10 option to ensure that
% only type 1 fonts will utilized at all point sizes. Without this option,
% it is possible that some math symbols, particularly those within
% footnotes, will be rendered in bitmap form which will result in a
% document that can not be IEEE Xplore compliant!
%
% Also, note that the amsmath package sets \interdisplaylinepenalty to 10000
% thus preventing page breaks from occurring within multiline equations. Use:
\interdisplaylinepenalty=2500
% after loading amsmath to restore such page breaks as IEEEtran.cls normally
% does. amsmath.sty is already installed on most LaTeX systems. The latest
% version and documentation can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/required/amslatex/math/





% *** SPECIALIZED LIST PACKAGES ***
%
%\usepackage{algorithmic}
% algorithmic.sty was written by Peter Williams and Rogerio Brito.
% This package provides an algorithmic environment fo describing algorithms.
% You can use the algorithmic environment in-text or within a figure
% environment to provide for a floating algorithm. Do NOT use the algorithm
% floating environment provided by algorithm.sty (by the same authors) or
% algorithm2e.sty (by Christophe Fiorio) as IEEE does not use dedicated
% algorithm float types and packages that provide these will not provide
% correct IEEE style captions. The latest version and documentation of
% algorithmic.sty can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/algorithms/
% There is also a support site at:
% http://algorithms.berlios.de/index.html
% Also of interest may be the (relatively newer and more customizable)
% algorithmicx.sty package by Szasz Janos:
% http://www.ctan.org/tex-archive/macros/latex/contrib/algorithmicx/




% *** ALIGNMENT PACKAGES ***
%
%\usepackage{array}
% Frank Mittelbach's and David Carlisle's array.sty patches and improves
% the standard LaTeX2e array and tabular environments to provide better
% appearance and additional user controls. As the default LaTeX2e table
% generation code is lacking to the point of almost being broken with
% respect to the quality of the end results, all users are strongly
% advised to use an enhanced (at the very least that provided by array.sty)
% set of table tools. array.sty is already installed on most systems. The
% latest version and documentation can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/required/tools/


%\usepackage{mdwmath}
%\usepackage{mdwtab}
% Also highly recommended is Mark Wooding's extremely powerful MDW tools,
% especially mdwmath.sty and mdwtab.sty which are used to format equations
% and tables, respectively. The MDWtools set is already installed on most
% LaTeX systems. The lastest version and documentation is available at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/mdwtools/


% IEEEtran contains the IEEEeqnarray family of commands that can be used to
% generate multiline equations as well as matrices, tables, etc., of high
% quality.


%\usepackage{eqparbox}
% Also of notable interest is Scott Pakin's eqparbox package for creating
% (automatically sized) equal width boxes - aka "natural width parboxes".
% Available at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/eqparbox/





% *** SUBFIGURE PACKAGES ***
%\usepackage[tight,footnotesize]{subfigure}
% subfigure.sty was written by Steven Douglas Cochran. This package makes it
% easy to put subfigures in your figures. e.g., "Figure 1a and 1b". For IEEE
% work, it is a good idea to load it with the tight package option to reduce
% the amount of white space around the subfigures. subfigure.sty is already
% installed on most LaTeX systems. The latest version and documentation can
% be obtained at:
% http://www.ctan.org/tex-archive/obsolete/macros/latex/contrib/subfigure/
% subfigure.sty has been superceeded by subfig.sty.



%\usepackage[caption=false]{caption}
%\usepackage[font=footnotesize]{subfig}
% subfig.sty, also written by Steven Douglas Cochran, is the modern
% replacement for subfigure.sty. However, subfig.sty requires and
% automatically loads Axel Sommerfeldt's caption.sty which will override
% IEEEtran.cls handling of captions and this will result in nonIEEE style
% figure/table captions. To prevent this problem, be sure and preload
% caption.sty with its "caption=false" package option. This is will preserve
% IEEEtran.cls handing of captions. Version 1.3 (2005/06/28) and later
% (recommended due to many improvements over 1.2) of subfig.sty supports
% the caption=false option directly:
\usepackage[caption=false,font=footnotesize]{subfig}
%
% The latest version and documentation can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/subfig/
% The latest version and documentation of caption.sty can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/caption/




% *** FLOAT PACKAGES ***
%
%\usepackage{fixltx2e}
% fixltx2e, the successor to the earlier fix2col.sty, was written by
% Frank Mittelbach and David Carlisle. This package corrects a few problems
% in the LaTeX2e kernel, the most notable of which is that in current
% LaTeX2e releases, the ordering of single and double column floats is not
% guaranteed to be preserved. Thus, an unpatched LaTeX2e can allow a
% single column figure to be placed prior to an earlier double column
% figure. The latest version and documentation can be found at:
% http://www.ctan.org/tex-archive/macros/latex/base/



%\usepackage{stfloats}
% stfloats.sty was written by Sigitas Tolusis. This package gives LaTeX2e
% the ability to do double column floats at the bottom of the page as well
% as the top. (e.g., "\begin{figure*}[!b]" is not normally possible in
% LaTeX2e). It also provides a command:
%\fnbelowfloat
% to enable the placement of footnotes below bottom floats (the standard
% LaTeX2e kernel puts them above bottom floats). This is an invasive package
% which rewrites many portions of the LaTeX2e float routines. It may not work
% with other packages that modify the LaTeX2e float routines. The latest
% version and documentation can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/sttools/
% Documentation is contained in the stfloats.sty comments as well as in the
% presfull.pdf file. Do not use the stfloats baselinefloat ability as IEEE
% does not allow \baselineskip to stretch. Authors submitting work to the
% IEEE should note that IEEE rarely uses double column equations and
% that authors should try to avoid such use. Do not be tempted to use the
% cuted.sty or midfloat.sty packages (also by Sigitas Tolusis) as IEEE does
% not format its papers in such ways.


%\ifCLASSOPTIONcaptionsoff
%  \usepackage[nomarkers]{endfloat}
% \let\MYoriglatexcaption\caption
% \renewcommand{\caption}[2][\relax]{\MYoriglatexcaption[#2]{#2}}
%\fi
% endfloat.sty was written by James Darrell McCauley and Jeff Goldberg.
% This package may be useful when used in conjunction with IEEEtran.cls'
% captionsoff option. Some IEEE journals/societies require that submissions
% have lists of figures/tables at the end of the paper and that
% figures/tables without any captions are placed on a page by themselves at
% the end of the document. If needed, the draftcls IEEEtran class option or
% \CLASSINPUTbaselinestretch interface can be used to increase the line
% spacing as well. Be sure and use the nomarkers option of endfloat to
% prevent endfloat from "marking" where the figures would have been placed
% in the text. The two hack lines of code above are a slight modification of
% that suggested by in the endfloat docs (section 8.3.1) to ensure that
% the full captions always appear in the list of figures/tables - even if
% the user used the short optional argument of \caption[]{}.
% IEEE papers do not typically make use of \caption[]'s optional argument,
% so this should not be an issue. A similar trick can be used to disable
% captions of packages such as subfig.sty that lack options to turn off
% the subcaptions:
% For subfig.sty:
% \let\MYorigsubfloat\subfloat
% \renewcommand{\subfloat}[2][\relax]{\MYorigsubfloat[]{#2}}
% For subfigure.sty:
% \let\MYorigsubfigure\subfigure
% \renewcommand{\subfigure}[2][\relax]{\MYorigsubfigure[]{#2}}
% However, the above trick will not work if both optional arguments of
% the \subfloat/subfig command are used. Furthermore, there needs to be a
% description of each subfigure *somewhere* and endfloat does not add
% subfigure captions to its list of figures. Thus, the best approach is to
% avoid the use of subfigure captions (many IEEE journals avoid them anyway)
% and instead reference/explain all the subfigures within the main caption.
% The latest version of endfloat.sty and its documentation can obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/endfloat/
%
% The IEEEtran \ifCLASSOPTIONcaptionsoff conditional can also be used
% later in the document, say, to conditionally put the References on a
% page by themselves.





% *** PDF, URL AND HYPERLINK PACKAGES ***
%
\usepackage{url}
% url.sty was written by Donald Arseneau. It provides better support for
% handling and breaking URLs. url.sty is already installed on most LaTeX
% systems. The latest version can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/misc/
% Read the url.sty source comments for usage information. Basically,
% \url{my_url_here}.





% *** Do not adjust lengths that control margins, column widths, etc. ***
% *** Do not use packages that alter fonts (such as pslatex).         ***
% There should be no need to do such things with IEEEtran.cls V1.6 and later.
% (Unless specifically asked to do so by the journal or conference you plan
% to submit to, of course. )
\usepackage[latin1]{inputenc} % f�r ���..
%\usepackage{amssymb}
\usepackage[squaren]{SIunits}
\usepackage{bm}
\usepackage{multirow}
\usepackage{color}
\usepackage[usenames,dvipsnames]{xcolor}
\usepackage{algorithm}
\usepackage{algorithmic}
\usepackage{abbrevs}
\usepackage{paralist}

% correct bad hyphenation here
\hyphenation{op-tical net-works semi-conduc-tor}

\newcommand{\notationObj}[1]{\hspace{-2mm}\begin{minipage}{0.9\columnwidth} \setlength{\parindent}{-15pt} $\bullet$ #1 \end{minipage}}

\newcommand{\comment}[1]{{\color{BrickRed}#1}}

\begin{document}
\setlength{\interdisplaylinepenalty}{2500}
% paper title
% can use linebreaks \\ within to get better formatting as desired
\title{An extended target tracking model with multiple random matrices and unified kinematics}
%
%
% author names and IEEE memberships
% note positions of commas and nonbreaking spaces ( ~ ) LaTeX will not break
% a structure at a ~ so this keeps an author's name from being broken across
% two lines.
% use \thanks{} to gain access to the first footnote area
% a separate \thanks must be used for each paragraph as LaTeX2e's \thanks
% was not built to handle multiple paragraphs
%

\author{Karl~Granstr\"om% <-this % stops a space ,~\IEEEmembership{Member,~IEEE}
%\thanks{M. Shell is with the Department of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, GA, 30332 USA e-mail: (see http://www.michaelshell.org/contact.html).}% <-this % stops a space
\thanks{Karl~Granstr\"om is with the Department
of Electrical Engineering, Division of Automatic Control, Link\"oping University, Link\"oping, SE-581~83, Sweden, e-mail: \texttt{karl@isy.liu.se}.}% <-this % stops a space
%\thanks{Manuscript received April 19, 2005; revised January 11, 2007.}}
}% note the % following the last \IEEEmembership and also \thanks -
% these prevent an unwanted space from occurring between the last author name
% and the end of the author line. i.e., if you had this:
%
% \author{....lastname \thanks{...} \thanks{...} }
%                     ^------------^------------^----Do not want these spaces!
%
% a space would be appended to the last name and could cause every name on that
% line to be shifted left slightly. This is one of those "LaTeX things". For
% instance, "\textbf{A} \textbf{B}" will typeset as "A B" not "AB". To get
% "AB" then you have to do: "\textbf{A}\textbf{B}"
% \thanks is no different in this regard, so shield the last } of each \thanks
% that ends a line with a % and do not let a space in before the next \thanks.
% Spaces after \IEEEmembership other than the last one are OK (and needed) as
% you are supposed to have spaces between the names. For what it is worth,
% this is a minor point as most people would not even notice if the said evil
% space somehow managed to creep in.



% The paper headers
%\markboth{Journal of \LaTeX\ Class Files,~Vol.~6, No.~1, January~2007}%
%{Shell \MakeLowercase{\textit{et al.}}: Bare Demo of IEEEtran.cls for Journals}
% The only time the second header will appear is for the odd numbered pages
% after the title page when using the twoside option.
%
% *** Note that you probably will NOT want to include the author's ***
% *** name in the headers of peer review papers.                   ***
% You can use \ifCLASSOPTIONpeerreview for conditional compilation here if
% you desire.




% If you want to put a publisher's ID mark on the page you can do it like
% this:
%\IEEEpubid{0000--0000/00\$00.00~\copyright~2007 IEEE}
% Remember, if you use this you must call \IEEEpubidadjcol in the second
% column for its text to clear the IEEEpubid mark.



% use for special paper notices
%\IEEEspecialpapernotice{(Invited Paper)}




% make the title area
\maketitle


\begin{abstract}
%\boldmath

This paper presents a model for tracking of extended targets, where each target is represented by a given number of elliptic subobjects. A gamma Gaussian inverse Wishart implementation is derived, and necessary approximations are suggested to alleviate the data association complexity. A simulation study shows the merits of the model compared to previous work on the topic.

\end{abstract}
% IEEEtran.cls defaults to using nonbold math in the Abstract.
% This preserves the distinction between vectors and scalars. However,
% if the journal you are submitting to favors bold math in the abstract,
% then you can use LaTeX's standard command \boldmath at the very start
% of the abstract to achieve this. Many IEEE journals frown on math
% in the abstract anyway.

% Note that keywords are not normally used for peerreview papers.
\begin{IEEEkeywords}
Target tracking, extended target, group target, measurement rate, random matrix, gamma distribution, Gaussian distribution, inverse Wishart distribution.
\end{IEEEkeywords}






% For peer review papers, you can put extra information on the cover
% page as needed:
% \ifCLASSOPTIONpeerreview
% \begin{center} \bfseries EDICS Category: 3-BBND \end{center}
% \fi
%
% For peerreview papers, this IEEEtran command inserts a page break and
% creates the second title. It will be ignored for other modes.
\IEEEpeerreviewmaketitle



%\section{Introduction}
% The very first letter is a 2 line initial drop letter followed
% by the rest of the first word in caps.
%
% form to use if the first word consists of a single letter:
% \IEEEPARstart{A}{demo} file is ....
%
% form to use if you need the single drop letter followed by
% normal text (unknown if ever used by IEEE):
% \IEEEPARstart{A}{}demo file is ....
%
% Some journals put the first two words in caps:
% \IEEEPARstart{T}{his demo} file is ....
%
% Here we have the typical use of a "T" for an initial drop letter
% and "HIS" in caps to complete the first word.
%\IEEEPARstart{T}{his} demo file is intended to serve as a ``starter file''
%for IEEE journal papers produced under \LaTeX\ using
%IEEEtran.cls version 1.7 and later.
% You must have at least 2 lines in the paragraph with the drop letter
% (should never be an issue)
%\IEEEPARstart{I}{n}

%\tableofcontents

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%
%%%%% Introduction
%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}
Target tracking can be defined as the processing of a sequence of measurements obtained from a target in order to maintain an estimate of the target's current state. In this context a point target is defined as a target which is assumed to give rise to at most one measurement per time step. With modern and more accurate sensors the target may occupy multiple resolution cells of the sensor, thus potentially giving rise to more than one measurement per time step. An extended target is defined as a target that potentially gives rise to more than one measurement per time step. Examples of extended target tracking include vehicle tracking using automotive radar and pedestrian tracking using laser range sensors. Closely related to extended target is group target, defined as a cluster of point targets which cannot be tracked individually, but has to be treated as a single object.

In point target tracking the estimated state typically corresponds to the targets position and its kinematics (velocity, heading, etc). In extended target tracking the multiple measurements make it possible to estimate also the target's extension in the measurement domain, \iep to estimate the shape, the size and the orientation of the target. To estimate the target's extension requires a measurement model that relates the multiple measurements to the states that govern the extension. 

Spatial distribution models in extended target tracking appeared in \cite{GilholmGMS:2005,GilholmS:2005}. Under this model each extended target measurement is a random sample from a probability distribution that is dependent on the extended target state. A number of different extended target models have been presented, where the targets are modeled as, \egp, sticks \cite{GilholmS:2005,BaumFH:2012,Boersetal:06a}, circles \cite{PetrovMGA:2011}, ellipses \cite{Koch:2008,BaumNH:2010,GranstromLO:2011,ZhuHL:2011,ReuterD:2011,DegermanWS:2011,LanRL:2012,ReuterWD:2012}, rectangles \cite{GranstromLO:2011}, or general shapes \cite{LundquistGO:11,BaumH:2011,LanRL:2014,LanRL:2012irregular}.

In this paper we consider state estimation for extended targets whose extensions cannot be approximated by a simple geometric shape such as an ellipse or a rectangle.
%Lan and Rong Li \cite{LanRL:2012irregular} give the following definition of an extended target with an irregular extension:
%\begin{definition}[Irregular Extended Target]
%An extended target whose extension cannot be simplified as an ellipsoid without losing important and useful extension information (\egp shape and orientation).\hfill$\square$
%\end{definition}
The extended target is modeled as a collection of elliptical subobjects, see \figurename~\ref{fig:multiple_ellipse_example_target}, and the positions and extensions of the subobjects are Gaussian inverse Wishart distributed. The scope of the paper is limited by the assumptions that
\begin{inparaenum}[\itshape a\upshape)]
	\item there is exactly one target present;
	\item there are no clutter measurements; and
	\item the number of subobjects is constant and known.
\end{inparaenum}
To handle multiple targets and clutter, the presented work can be integrated into a multiple target framework, \egp an extended target \phd/\cphd filter \cite{mahler_FUSION_2009_extTarg,GranstromLO:2010,GranstromLO:2012,GranstromO:2012a,LundquistGO:2013,SwainC:2010,SwainC:2012}. Estimating the number of subobjects is left for future work.

%The rest of the paper is organized as follows. In the next section we give an overview of the random matrix framework for extended target modeling, and list the contributions of the paper. In Section~\ref{sec:new_multi_ellipse_model} the extended target model is presented, and in Section~\ref{sec:GGIW_implementation} a gamma Gaussian inverse Wishart implementation is given. Metrics for performance evaluation are given in Section~\ref{sec:performance_eval}, and results from a simulation study are presented in Section~\ref{sec:simulation_results}. The paper is concluded in Section~\ref{sec:conclusion_future_work}.

\begin{figure}[t]
\centering
\subfloat{\includegraphics[width=0.40\columnwidth]{multiple_ellipse_example_target}}%
\hfil
\subfloat{\includegraphics[width=0.40\columnwidth]{multiple_ellipse_example_target_Vshape}}%
\caption{Examples in 2D of extended/group targets that are represented by elliptic subobjects. Neither one of the examples has a shape that can be described by a single ellipse. Left: The overall position of the extended target is $\mathbf{p}$, and coincides with one of the subobject's position. The positions of remaining subobjects are given by offsets $\mathbf{d}$ from the overall position. Right: The overall position does not coincide with either subobject's position.}
\label{fig:multiple_ellipse_example_target}
\end{figure}

\begin{table}
\caption{Notation}
\label{tab:notation}
\vspace{-5mm}
\rule[0pt]{\columnwidth}{1pt}
$\bullet \ $ $\mathbb{R}^{n}$ is the set of real column vectors of length $n$, $\mathbb{S}_{++}^{n}$ is the set of symmetric positive definite $n\times n$ matrices, $\mathbb{S}_{+}^{n}$ is the set of symmetric positive semi-definite $n\times n$ matrices, and $\mathbb{N}$ is the set of non-negative integers.	

$\bullet \ $ $\Id$ is a $d\times d$ eye matrix, and $\mathbf{0}_{d\times e}$ is a $d\times e$ all-zero matrix.

$\bullet \ $ $|\cdot|$ is absolute value, $\left\|\cdot\right\|_{2}$ is Euclidean norm, and $\left\|\cdot\right\|_{F}$ is Frobenius norm.

$\bullet \ $ $\mathcal{PS}\left(n;\ \gamma\right)$ denotes a Poisson probability mass function (pmf) defined of the integer $n\in\mathbb{N}$ with rate parameter $\gamma>0$,
		\begin{align*}
			\mathcal{PS}\left(n;\ \gamma\right) = & {\gamma^{n}e^{-\gamma}}{\left(n!\right)^{-1}}.
		\end{align*}	

$\bullet \ $ $\Gammapdf{\gamma}{\alpha}{\beta}$ denotes a gamma probability density function (pdf) defined over the scalar $\gamma >0$ with scalar shape parameter $\alpha > 0$ and scalar inverse scale parameter $\beta > 0$,
\begin{align*}
	\Gammapdf{\gamma}{\alpha}{\beta} = {\beta^{\alpha}}{\Gamma(\alpha)^{-1}}\gamma^{\alpha-1}e^{-\beta\gamma},
\end{align*}
where $\Gamma(\cdot)$ is the gamma function.

$\bullet \ $ $\Npdfbig{\sx}{\GaussMean}{P}$ denotes a multi-variate Gaussian pdf defined over the vector $\sx\in\mathbb{R}^{n_x}$ with mean vector $\GaussMean\in\mathbb{R}^{n_x}$, and covariance matrix $P\in\mathbb{S}_{+}^{n_x}$,
	\begin{align*}
	\Npdfbig{\sx}{\GaussMean}{P} = \frac{\exp\left({-\frac{1}{2}\left(\sx-\GaussMean\right)^{\tp}P^{-1}\left(\sx-\GaussMean\right)}\right)}{\left(2\pi\right)^{\frac{n_x}{2}}\det(P)^{\frac{1}{2}}}.
	\end{align*}
	where $\det (\cdot)$ is the matrix determinant function.

$\bullet \ $ $\IWishpdf{\ext}{v}{V}$ denotes an inverse Wishart pdf defined over the matrix $\ext\in\mathbb{S}_{++}^{d}$ with scalar degrees of freedom $v>2d$ and parameter matrix $V\in\mathbb{S}_{++}^{d}$, \cite[Definition 3.4.1]{GuptaN:2000}
	\begin{align*}
	\IWishpdf{\ext}{v}{V} = \frac{2^{-\frac{v-d-1}{2}}\det(V)^{\frac{v-d-1}{2}}}{\Gamma_{d}\left(\frac{v-d-1}{2}\right)\det(\ext)^{\frac{v}{2}}}\etr\left(-\frac{1}{2}\ext^{-1}V\right),
	\end{align*}
	where $\etr(\cdot) = \exp \left(\tr(\cdot)\right)$ is exponential of the matrix trace, and $\Gamma_{d}\left(\cdot\right)$ is the multivariate gamma function. The multivariate gamma function can be expressed as a product of ordinary gamma functions, see \cite[Theorem 1.4.1]{GuptaN:2000}.

$\bullet \ $ $\Wishpdf{\ext}{w}{W}$ denotes a Wishart pdf defined over the matrix $\ext\in\mathbb{S}_{++}^{d}$ with scalar degrees of freedom $w\geq d$ and parameter matrix $W\in\mathbb{S}_{++}^{d}$, \cite[Definition 3.2.1]{GuptaN:2000}
	\begin{align*}
		\Wishpdf{\ext}{w}{W} = \frac{2^{-\frac{wd}{2}}\det(\ext)^{\frac{w-d-1}{2}}}{\Gamma_{d}\left(\frac{w}{2}\right)\det(W)^{\frac{n}{2}}}\etr\left(-\frac{1}{2}W^{-1}\ext\right).
	\end{align*}
\rule[0pt]{\columnwidth}{1pt}
\end{table}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%
%%%%% Overview of random matrix framework
%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Previous work and paper contributions}
%The presented extended target model belongs to a class of extended target models that can be called the random matrix framework. 
% In this section we will first review previous work in this framework, and we will then list the main contributions. %explain how the presented model relates to, and extends, previous work.

\subsection{Overview of random matrix framework}
Notation is given in Table~\ref{tab:notation}. In the random matrix extended target model, originally proposed by Koch in \cite{Koch:2008}, the extended target state is the combination of a kinematic state vector $\sx_{k}$ and an extension matrix $\ext_{k}$. The vector $\sx_{k}$ represents the target's position and kinematics, and the matrix $\ext_{k}$ represents the target's size and shape, \iep its spatial extension. The matrix $\ext_{k}$ is modeled as being symmetric and positive definite, which implies that the target shape is approximated by an ellipse. The ellipse shape may seem limiting, however the model is applicable to many real scenarios, e.g. pedestrian tracking \cite{GranstromO:2012a}.

In \cite{Koch:2008} the target state, the target generated measurements, and the transition density, are modeled as
\begin{subequations}
\begin{align}
p\left(\sx_{k},\ext_{k}|\setZ^{k}\right) = & p\left(\sx_{k}|\ext_{k},\setZ^{k}\right)p\left(\ext_{k}|\setZ^{k}\right) \\
= & \Npdfbig{\sx_{k}}{m_{k|k}}{P_{k|k}\otimes\ext_{k}} \nonumber \\
& \times \IWishpdf{\ext_{k}}{v_{k|k}}{V_{k|k}}, \label{eq:Koch_GIW_model}\\
p\left(\sz_{k}|\sx_{k},\ext_{k}\right) = & \Npdfbig{\sz_{k}}{H_{k}\sx_{k}}{\ext_{k}}. \label{eq:Linear_Gaussian_Measurements}\\
p\left(\sx_{k+1},\ext_{k+1}|\sx_{k},\ext_{k}\right) = & p\left(\sx_{k+1}|\ext_{k+1},\sx_{k}\right)p\left(\ext_{k+1}|\ext_{k}\right), \label{eq:Koch_transition_density}
\end{align}%
\label{eq:Koch_random_matrix_model}%
\end{subequations}%
In this model the kinematic state $\sx_{k}$ consists of a spatial state component $\mathbf{r}_{k}$  (target position), and derivatives of $\mathbf{r}_{k}$ (typically velocity and acceleration) \cite{Koch:2008}. Non-linear dynamics, such as turn-rate, are not included in the kinematic vector. The measurement update is linear \cite{Koch:2008}, a derivation of the predicted likelihood can be found in \cite[Appendix A]{GranstromO:2012a}. A linear Gaussian transition density is used for the kinematic state, and for the extension a simple heuristic is used in which the expected value is kept constant and the variance is increased \cite{Koch:2008}. The extension transition density $p\left(\ext_{k+1}|\ext_{k}\right)$ in \eqref{eq:Koch_transition_density} assumes independence of the kinematic state $\sx_{k}$, which does not account for, e.g., rotations during a turning maneuver \cite{GranstromO:2014}.

The random matrix model \eqref{eq:Koch_random_matrix_model} was modified in \cite{FeldmannF:2008,FeldmannFK:2011}, where the target state, the target generated measurements, and the transition density, are modeled as
\begin{subequations}
\begin{align}
p\left(\sx_{k},\ext_{k}|\setZ^{k}\right)  \approx &  p\left(\sx_{k}|\setZ^{k}\right)p\left(\ext_{k}|\setZ^{k}\right) \label{eq:Feldmann_random_matrix_state} \\
= & \Npdfbig{\sx_{k}}{m_{k|k}}{P_{k|k}} \nonumber \\
& \times\IWishpdf{\ext_{k}}{v_{k|k}}{V_{k|k}},\\
p\left(\sz_{k}|\sx_{k},\ext_{k}\right) = & \Npdfbig{\sz_{k}}{H_{k}\sx_{k}}{z\ext_{k}+R}, \label{eq:Feldmann_measurement_model} \\
p\left(\sx_{k+1},\ext_{k+1}|\sx_{k},\ext_{k}\right) = & p\left(\sx_{k+1}|\sx_{k}\right)p\left(\ext_{k+1}|\ext_{k}\right), \label{eq:Feldmann_transition_density}
\end{align}%
\label{eq:Feldmann_random_matrix_model}%
\end{subequations}%
where $z$ is a scaling factor and $R$ is measurement noise. Note the assumed independence between the kinematic state $\sx_{k}$ and $\ext_{k}$ in \eqref{eq:Feldmann_random_matrix_state}, an assumption that cannot be fully theoretically justified\footnote{Conditioned on a set of measurements $\setZ$ the kinematic state $\sx$ and extension state $\ext$ are necessarily dependent.}. Further the measurement update is no longer linear and must be approximated, see \cite{FeldmannFK:2011} for details. However, there are considerable practical advantages to the model \cite{FeldmannFK:2011}.

This model allows for a more general class of kinematic state vectors $\sx_{k}$, e.g. including non-linear dynamics such as heading and turn-rate, and the Gaussian covariance is no longer intertwined with the extension matrix. This measurement model is better when the size of the extension and the size of the sensor noise are within the same order of magnitude \cite{FeldmannFK:2011}. The assumed independence between $\sx_{k}$ and $\ext_{k}$ is alleviated in practice by the measurement update which provides for the necessary interdependence between kinematics and extension estimation \cite{FeldmannFK:2011}. An alternative measurement update for the measurement model \eqref{eq:Feldmann_measurement_model}, based on variational Bayes approximation, is given in \cite{Orguner:2012}.

The kinematics transition density $p\left(\sx_{k+1}|\sx_{k}\right)$ in \eqref{eq:Feldmann_transition_density} is assumed independent of the extension. This neglects factors such as wind resistance, which can be modeled as a function of the extension $\ext_{k}$, however the assumption is necessary to retain the functional form \eqref{eq:Feldmann_random_matrix_state} in a Bayesian recursion. A linear Gaussian transition density is used for the kinematic state, and a heuristic transition similar to the one in \cite{Koch:2008} is used for the extension.

An alternative to the heuristic extension predictions from \cite{Koch:2008,FeldmannFK:2011} is to use a Wishart transition density \cite{Koch:2008}, see also \cite{LianHLYZ:2010,LanRL:2012,GranstromO:2014}. In \cite{LanRL:2012} transformations of the extension are allowed via parameter matrices $A_{k}$,
\begin{align}
p(\ext_{k+1}|\ext_{k})=&\Wishpdf{\ext_{k+1}}{\delta_{k}}{A_{k}\ext_{k}A_{k}^{\tp}}. \label{eq:LanRongLi_prediction}
\end{align}
The parameter matrices correspond to, \egp, rotation matrices. This is generalized in \cite{GranstromO:2014} to allow for transformation matrices $M(\sx_{k})$ that are functions of the kinematic state,
\begin{align}
p(\ext_{k+1}|\sx_{k},\ext_{k}) = & \Wishpdf{\ext_{k+1}}{n_{k}}{\frac{M({\sx_{k}})\ext_{k}(M({\sx_{k}}))^{\tp}}{n_{k}}}, \label{eq:GranstromOrguner_prediction}
\end{align}
which means that the rotation angle can be coupled to, e.g., the turn-rate and estimated online. The transition density \eqref{eq:GranstromOrguner_prediction} relaxes the assumption (made in \eqref{eq:Koch_transition_density}, \eqref{eq:Feldmann_transition_density}, and \eqref{eq:LanRongLi_prediction}) that the extension's time evolution is independent of the kinematic state. A comparison of the models \eqref{eq:Koch_transition_density}, \eqref{eq:Feldmann_transition_density}, \eqref{eq:LanRongLi_prediction} and \eqref{eq:GranstromOrguner_prediction} is given in \cite{GranstromO:2014}, where \eqref{eq:GranstromOrguner_prediction} is shown to give lower errors at lower computational complexity.

In addition to the transition density \eqref{eq:LanRongLi_prediction}, a measurement model is also suggested in \cite{LanRL:2012},
\begin{align}
p\left(\sz_{k}|\sx_{k},\ext_{k}\right) = & \Npdfbig{\sz_{k}}{H_{k}\sx_{k}}{B_{k}\ext_{k}B_{k}^{\tp}} \label{eq:LanRL_measurement_model}
\end{align}
where $B_{k}$ is a parameter matrix. Under the assumption $\ext_{k}\approx\hat{\ext}_{k|k-1}=\E[\ext_{k}|\setZ^{k-1}]$  the model \eqref{eq:LanRL_measurement_model} incorporates \eqref{eq:Feldmann_measurement_model} approximately when $B_{k} = (z\hat{\ext}_{k|k-1}+R)^{1/2}\hat{\ext}_{k|k-1}^{-1/2}$.

The random matrix model has been integrated into the Probabilistic Multi-Hypothesis Tracking (\textsc{pmht}) framework \cite{StreitL:1993}, see \cite{WienekeK:2010,WienekeD:2011,WienekeK:2012}. The model has also been used in \phd- and \cphd-filters for multiple extended target tracking in clutter, see \cite{GranstromO:2012a,LundquistGO:2013}.

In \cite{LanRL:2014} a single extended target model is given where the extended target is a combination of multiple subobjects with kinematic state vectors $\sx_{k}^{(i)}$ and extension matrices $\ext_{k}^{(i)}$. Each subobject is modeled using \eqref{eq:Koch_GIW_model},
\begin{align}
\Npdfbig{\sx_{k}^{(i)}}{m_{k|k}^{(i)}}{P_{k|k}^{(i)}\otimes\ext_{k}^{(i)}}\IWishpdf{\ext_{k}^{(i)}}{v_{k|k}^{(i)}}{V_{k|k}^{(i)}}. \label{eq:LanRongLiStateModel}
\end{align}
%and the measurement model \eqref{eq:LanRL_measurement_model} and transition density \eqref{eq:LanRongLi_prediction} are used. 
Using multiple instances of a simpler shape alleviates the limitations posed by the implied elliptic target shape\footnote{As the number of ellipses grows, their combination can form nearly any given shape.}, and also retains, on a subobject level, the simplicity of the random matrix model \cite{Koch:2008,FeldmannFK:2011}.

% In \cite{LanRL:2014} a multiple model approach is used to model different types of motion, where the subobject motion is split into both individual and common parts. Using the Wishart transition density \eqref{eq:LanRongLi_prediction} common kinematics for all extension matrices $\ext_{k}^{(i)}$ is achieved using the same parameter matrix for all subobjects. However, there are no common kinematics for the subobject kinematic state vectors $\sx_{k}^{(i)}$.

\subsection{Paper contributions}
The major contributions in this paper are:

$\bullet$ A new state space model where all subobjects's positions are modeled as fully dependent using a single state vector, there are fully unified kinematics, and where the measurement rate and extension of each subobject are individually modeled.

$\bullet$ A derivation of the prediction update, the measurement update, and the predicted likelihood.

$\bullet$ A computationally effcient gamma Gaussian inverse Wishart implementation, including an initialization method that does not rely on any a priori information about the target.

$\bullet$ A simple and effective method that minimizes the number of association events that have to be considered without relying on a predicted target estimate. The method is capable of handling partial occlusion of the extended target, i.e. one or a few of the subobjects are hidden from sensor view.

$\bullet$ The proposed extended target model is validated on simulated data from realistic scenarios, and the results are compared to previous work on the topic. %Additionally, if the paper is accepted for publication, Matlab code will be made available via the author's homepage.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%
%%%%% Proposed multi ellipse model
%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Proposed multiple ellipse model}
\label{sec:new_multi_ellipse_model}

\subsection{Extended target state}
The extended target is made up of a combination of $N_{s,k}$ $d$-dimensional subobjects, where $N_{s,k}$ is known. Each subobject $i$ is described by a position $\mathbf{p}_{k}^{(i)}\in\mathbb{R}^{d}$, a measurement rate $\gamma_{k}^{(i)}>0$ and an extension state $\ext_{k}^{(i)}\in\mathbb{S}_{++}^{d}$, where sub-index $k$ refers to discrete time step $t_k$. The measurement rate governs how many measurements the subobject generates per time step, and the extension describes the size and the shape of the subobject. Because extended targets in most cases can be assumed to be rigid bodies the subobjects have unified dynamics, by which we mean that all subobjects move forward with the same velocity and the same heading, turn with the same turn-rate, etc. The unified dynamics are denoted $\mathbf{c}_{k}\in\mathbb{R}^{n_c}$, where $\mathbf{c}_{k}$ includes parameters for, \egp, velocity, acceleration, heading and turn-rate. Note that $\mathbf{c}_{k}$, in addition to parameters for unified dynamics, also may include parameters for individual subobject dynamics. This is useful for group tracking, where the individual targets in the group may shift their positions within the group.

The positions $\mathbf{p}_{k}^{(i)}$ of the subobjects are $\mathbf{p}_{k}^{(1)} = \mathbf{p}_{k}$ and $\mathbf{p}_{k}^{(i)} = \mathbf{p}_{k} + \mathbf{d}_{k}^{(i)}$, $i=2,\ldots,N_{s,k}$,
%\begin{subequations}
%\begin{align}
%\mathbf{p}_{k}^{(1)} & = \mathbf{p}_{k}, \\
%\mathbf{p}_{k}^{(i)} & = \mathbf{p}_{k} + \mathbf{d}_{k}^{(i)}, \quad i=2,\ldots,N_{s,k},
%\end{align}
%\end{subequations}
\iep the positions of subobjects $i=2,\ldots,N_{s,k}$ are offset by vectors $\mathbf{d}_{k}^{(i)}\in\mathbb{R}^{d}$ from the first subobject's position $\mathbf{p}_{k}\in\mathbb{R}^{d}$. The first subobject is also referred to as the main subobject, and the position $\mathbf{p}_{k}$ is referred to as the overall position. The unified kinematics are defined w.r.t. the overall position. For linear dynamics it does not matter which subobject is denoted the first subobject. However, this is important for non-linear dynamics, \egp a turning maneuver that causes the extended target to rotate.

The positions and dynamics of all subobjects are jointly described by a kinematic state $\sx_{k}\in\mathbb{R}^{n_x}$,
\begin{align}
\sx_{k} = \begin{bmatrix} \mathbf{p}_{k}^{\tp} & \mathbf{c}_{k}^{\tp} & \left(\mathbf{d}_{k}^{(2)}\right)^{\tp} & \ldots &  \left(\mathbf{d}_{k}^{(N_{s,k})}\right)^{\tp} \end{bmatrix}^{\tp}. \label{eq:kinematic_state_vector}
\end{align}
Note that the position of one of the subobjects must coincide with the overall position, because if $\mathbf{p}_{k}^{(i)} = \mathbf{p}_{k} + \mathbf{d}_{k}^{(i)}$ for all $i$ then $\mathbf{p}_{k}$ is not observable.

For brevity the measurement rates, kinematic state and extension states are abbreviated as follows
\begin{subequations}
\begin{align}
\xi_{k} = & \left(\gamma_{k}^{(1)},\ldots,\gamma_{k}^{(N_{s,k})},\sx_{k},\ext_{k}^{(1)},\ldots,\ext_{k}^{(N_{s,k})}\right) \\
% = & \left(\left\{\gamma_{k}^{(i)}\right\}_{i=1}^{N_{s,k}},\sx_{k},\left\{\ext_{k}^{(i)}\right\}_{i=1}^{N_{s,k}}\right) \\
= & \left(\setGamma_{k},\sx_{k},\setX_{k}\right)
\end{align}%
\end{subequations}%
where $\xi_{k}$ is referred to as the extended target state. Let $\setZ_{k}$ be a set of target generated measurements $\setZ_{k} = \{\sz_{k}^{(j)}\}_{j=1}^{n_{z,k}}$, $\sz_{k}^{(j)}\in\mathbb{R}^{d}, \ \forall j$,
%\begin{align}
%\setZ_{k} = \left\{\sz_{k}^{(j)}\right\}_{j=1}^{n_{z,k}}, \qquad \sz_{k}^{(j)}\in\mathbb{R}^{d}, \ \forall j,
%\end{align}
and let $\setZ^{k}$ be a sequence of measurement sets from time $t_{0}$ to time $t_{k}$. The distribution of the extended target state $\xi_{k}$, conditioned on the history of measurement sets, is
%\begin{subequations}
\begin{align}
%p\left(\xi_{k}\left|\setZ^{k}\right.\right) = & p\left(\setGamma_{k},\sx_{k},\setX_{k}\left|\setZ^{k}\right.\right) \\
p\left(\xi_{k}\left|\setZ^{k}\right.\right)= & p\left(\setGamma_{k}\left|\sx_{k},\setX_{k},\setZ^{k}\right.\right) p\left(\sx_{k}\left|\setX_{k},\setZ^{k}\right.\right) p\left(\setX_{k}\left|\setZ^{k}\right.\right),
\end{align}%
%\end{subequations}%
The following assumptions are made about the state $\xi_{k}$.
\begin{assumption}
The measurement rates are independent of the kinematic vector and the random matrices.\hfill$\square$
\end{assumption}
\begin{remark}
In reality the measurement rate is often dependent on the distance to the sensor (i.e. on the position) and on the size of the target (i.e. the extension). However, the distribution of the measurement rate, conditioned on the kinematic and extension states, is unknown in many applications. The variance of the estimated measurement rate is sufficient to model the variations over time \cite{GranstromO:2012c}, and the assumption ensures practical computational tractability. 
\hfill$\square$\end{remark}
\begin{assumption}
The kinematic vector is approximated as independent of the random matrices.\hfill$\square$
\end{assumption}
\begin{remark}
This assumption is analogous to \eqref{eq:Feldmann_random_matrix_state}, and similarly it neglects dependence between the subobjects positions and extension states -- an assumption that cannot be fully theoretically justified. However, the assumption is necessary to enable the subobject positions and unified kinematics to be modeled as a single state vector. Just as in \cite{FeldmannFK:2011,GranstromO:2014}, the time update and measurement update that are derived in this paper will provide for the necessary practical interdependency between the kinematic state estimate and the extension state estimate. 
\hfill$\square$\end{remark}
\begin{assumption}
The measurement rates are independent of each other. The same holds for the random matrices.\hfill$\square$
\end{assumption}
\begin{remark}
For symmetric targets, \egp airplanes, the subobjects that correspond to the wings are not fully independent since the wings are symmetric. Modeling how random matrices, or measurement rates, depend on each other is difficult in a general case, and the assumption simplifies further analysis greatly and ensures practical computational tractability.  An important topic for future work is to consider how symmetry can be used to relax this assumption. For the random matrices, the same is assumed in \cite{LanRL:2014}, see \eqref{eq:LanRongLiStateModel}. 
\hfill$\square$\end{remark}

This gives the following extended target state distribution,
%\begin{subequations}
\begin{align}
%& p \left(\xi_{k} \left| \setZ^{k}\right.\right) \nonumber \\
% = & p\left(\setGamma_{k}\left|\mode_{k},\setZ^{k}\right.\right) p\left(\sx_{k}\left|\mode_{k},\setZ^{k}\right.\right) p\left(\setX_{k}\left|\mode_{k},\setZ^{k}\right.\right) p\left(\mode_{k}\left|\setZ^{k}\right.\right) \\
p \left(\xi_{k} \left| \setZ^{k}\right.\right) = & p\left(\sx_{k}\left|\setZ^{k}\right.\right) \prod_{i=1}^{N_{s,k}} p\left(\gamma_{k}^{(i)} \left| \setZ^{k}\right.\right) p\left(\ext_{k}^{(i)} \left| \setZ^{k} \right.\right) . \label{eq:multiple_ellipse_ET_distribution}%
\end{align}%
%\label{eq:multiple_ellipse_ET_distribution}%
%\end{subequations}%
%Note that while the dependence between the kinematic state and the random matrices is not explicitly modeled, the prediction and correction presented below provides for interdependency between the estimated quantities.

Because of the many uncertainties involved the extended target state distribution can also be represented by a distribution mixture. In this case the distribution is
\begin{align}
p \left(\xi_{k} \left| \setZ^{k}\right.\right) = \sum_{\ell=1}^{J} w^{(\ell)} p^{(\ell)} \left(\xi_{k} \left| \setZ^{k}\right.\right),
\end{align}
where $J$ is the number of components in the mixture, the weights $w^{(\ell)}$ sum to unity, and each $p^{(\ell)} \left(\cdot\right)$ are of the form \eqref{eq:multiple_ellipse_ET_distribution}. The weights can be interpreted as probabilities that the $\ell$:th mixture component $p^{(\ell)} \left(\cdot\right)$ is the true distribution.

\subsection{Prediction}
For brevity and increased readability, in this section we drop sub-index $_k$ and write sub-index $_{k+1}$ as sub-index $_{+}$, \iep we write $\xi$ and $\xi_{+}$ instead of $\xi_{k}$ and $\xi_{k+1}$. The state transition density $p\left(\xi_{+}|\xi_{}\right)$ describes the time evolution of the extended target state from time $t_{}$ to time $t_{+}$. The transition density decomposes as follows
%\begin{subequations}
\begin{align}
%p\left(\xi_{+}|\xi_{}\right) = & p\left( \setGamma_{+},\sx_{+},\setX_{+}\left|\setGamma_{},\sx_{},\setX_{} \right. \right) \\
%= & p\left( \setGamma_{+} \left|\sx_{+},\setX_{+},\setGamma_{},\sx_{},\setX_{} \right. \right) p\left(\sx_{+} \left|\setX_{+}, \setGamma_{},\sx_{},\setX_{} \right.\right) \nonumber \\
%& \times p\left(\setX_{+} \left|\setGamma_{},\sx_{},\setX_{} \right.\right) \\
p\left(\xi_{+}|\xi_{}\right) = & p\left( \setGamma_{+} \left|\sx_{+},\setX_{+},\setGamma_{} \right. \right) p\left(\sx_{+} \left|\setX_{+}, \setGamma_{},\sx_{} \right.\right) \nonumber \\
& \times p\left(\setX_{+} \left|\setGamma_{},\sx_{},\setX_{} \right.\right),
\end{align}
%\end{subequations}
where we have used Bayes rule and Markov-property assumptions. We now make the following assumptions:
\begin{assumption}
The measurement rates can be predicted independently of the kinematic vector and the random matrices. \hfill$\square$
\end{assumption}
\begin{remark}
This assumption is also made in \cite{GranstromO:2012c}. In reality the measurement rates typically depend on both the size of the target, and its distance from the sensor (\iep depends on the position), and this assumption neglects such dependencies. However, constructing a general model for these dependencies seems difficult, and the assumption simplifies further analysis significantly. Furthermore, the variances of the estimated measurement rates are sufficient to capture the variations over time \cite{GranstromO:2012c}.
\hfill$\square$\end{remark}
\begin{assumption}
The kinematic vector and the random matrices can be predicted independently of the prior measurement rates $\setGamma_{}$.\hfill$\square$
\end{assumption}
\begin{remark}
This assumption can be justified in most -- if not all -- practical cases, since neither the kinematics nor the size and shape of a target evolve differently depending on how many sensor detections the target generates.
\hfill$\square$\end{remark}
\begin{assumption}
The kinematic vector can be predicted independently of the random matrices.\hfill$\square$
\end{assumption}
\begin{remark}
This assumption is analogous to \eqref{eq:Feldmann_transition_density}, and similarly it neglects aspects such as wind resistance, which can be modeled as dependent on the target's size (\iep dependent on the random matrices). In this work the assumption is necessary for the predicted state distribution to be of the same functional form as the posterior state distribution, which is a typical requirement in Bayesian estimation. Both the measurement update and the prediction update that are used in this work provide for interdependency between the predicted kinematic vector and the random matrices. For single ellipse targets it is shown in \cite{GranstromO:2014} that estimation performance is not negatively affected by the assumption.
\hfill$\square$\end{remark}
\begin{assumption}
Each measurement rate can be predicted independently of the other measurement rates. The same holds for the random matrices.\hfill$\square$
\end{assumption}
\begin{remark}
In practice the measurement rates and extensions evolve over time dependent on the position and kinematics. Because the rates and extensions are parts of the same object (i.e. subobjects), the evolution over time should be dependent. Under this assumption this dependence is not modeled, however the assumption is necessary for the predicted state distribution to be of the same functional form as the posterior. For the measurement rates the same motion model is used for all estimated rates, and the estimated variances are sufficient to model the time variations. For the extension states, the kinematic state is used in the prediction, which introduces sufficient dependence.
\hfill$\square$\end{remark}
%Following Assumptions~\ref{ass:MeasRate_pred_independently}, \ref{ass:KinVec_RM_pred_ind_MeasRate}, \ref{ass:KinVec_predInd_RM} and \ref{ass:MeasRate_RM_indPred} the transition density becomes
This gives the following transition density for $\xi_{}$,
%\begin{subequations}
\begin{align}
%p \left(\xi_{+}|\xi_{}\right) = & p\left( \setGamma_{+} \left| \setGamma_{} \right. \right) p\left(\sx_{+} \left| \sx_{} \right.\right) p\left(\setX_{+} \left| \sx_{},\setX_{} \right.\right) \\
p \left(\xi_{+}|\xi_{}\right) = & p\left(\sx_{+} \left| \sx_{} \right.\right) \prod_{i=1}^{N_{s,k}} p\left( \gamma_{+}^{(i)} \left| \gamma_{}^{(i)} \right. \right) p\left( \ext_{+}^{(i)} \left| \sx_{}, \ext_{}^{(i)} \right. \right) . \label{eq:multiple_ellipse_transition_density}
\end{align}%
%\label{eq:multiple_ellipse_transition_density}%
%\end{subequations}%
%In Bayesian state estimation the prediction step consists of solving the following equation,
%\begin{align}
%p\left(\xi_{+}\left|\setZ^{k}\right.\right) = & \int p(\xi_{+}|\xi_k) p\left(\xi_k|\setZ^{k}\right) \diff \xi_{}.
%\label{eq:prediction_integral}
%\end{align}
With a posterior \eqref{eq:multiple_ellipse_ET_distribution} and a transition density \eqref{eq:multiple_ellipse_transition_density} the Bayes predicted distribution is
\begin{subequations}
\begin{align}
p & \left(\xi_{+}\left|\setZ^{k}\right.\right) = \int p(\xi_{+}|\xi) p\left(\xi|\setZ^{k}\right) \diff \xi_{} \\
%\end{align}
%\begin{align}
%= & \iiint p\left( \setGamma_{+} \left| \setGamma_{} \right. \right) p\left(\sx_{+} \left| \sx_{} \right.\right) p\left(\setX_{+} \left| \sx_{},\setX_{} \right.\right) \nonumber \\
%& \times p\left(\setGamma_{}\left|\setZ^{k}\right.\right) p\left(\sx_{}\left|\setZ^{k}\right.\right) p\left(\setX_{}\left|\setZ^{k}\right.\right) \diff\setGamma_{}\diff\sx_{}\diff\setX_{} \\
%= & \int p\left( \setGamma_{+} \left| \setGamma_{} \right. \right) p\left(\setGamma_{}\left|\setZ^{k}\right.\right) \diff\setGamma_{} \\
%& \times \iint p\left(\sx_{+} \left| \sx_{} \right.\right) p\left(\setX_{+} \left| \sx_{},\setX_{} \right.\right) \nonumber \\
%& \times p\left(\sx_{}\left|\setZ^{k}\right.\right) p\left(\setX_{}\left|\setZ^{k}\right.\right) \diff\setX_{}\diff\sx_{} \nonumber \\
= & \left(\prod_{i=1}^{N_{s,k}}\int p\left( \gamma_{+}^{(i)} \left| \gamma_{}^{(i)} \right. \right) p\left(\gamma_{}^{(i)}\left|\setZ^{k}\right.\right) \diff\gamma_{}^{(i)}\right) \\
& \times \int \left( \prod_{i=1}^{N_{s,k}}\int p\left(\ext_{+}^{(i)} \left| \sx_{},\ext_{}^{(i)} \right.\right) p\left(\ext_{}^{(i)} \left| \setZ^{k}\right.\right) \diff\ext_{}^{(i)} \right) \nonumber \\
& \hspace{10mm}  \times p\left(\sx_{+} \left| \sx_{} \right.\right) p\left(\sx_{}\left|\setZ^{k}\right.\right) \diff\sx_{} \nonumber \\
= & \left(\prod_{i=1}^{N_{s,k}} p\left( \gamma_{+}^{(i)}\left|\setZ^{k}\right.\right) \right) \\
& \times \int \left( \prod_{i=1}^{N_{s,k}} p\left(\ext_{+}^{(i)} \left| \sx_{} , \setZ^{k} \right.\right) \right) p\left(\sx_{+} \left| \sx_{} \right.\right) p\left(\sx_{}\left|\setZ^{k}\right.\right)  \diff\sx_{} \nonumber
\end{align}
\end{subequations}
We want the predicted distribution to be of the same functional form as the posterior \eqref{eq:multiple_ellipse_ET_distribution}, however in general the following equality does not hold,
\begin{align}
\int & \left( \prod_{i=1}^{N_{s,k}} p\left(\ext_{+}^{(i)} \left| \sx_{} , \setZ^{k} \right.\right) \right) p\left(\sx_{+} \left| \sx_{} \right.\right) p\left(\sx_{}\left|\setZ^{k}\right.\right)  \diff\sx_{} \nonumber \\
= & p\left(\sx_{+}\left|\setZ^{k}\right.\right) \left( \prod_{i=1}^{N_{s,k}} p\left(\ext_{+}^{(i)} \left| \setZ^{k} \right.\right) \right).
\end{align}
Therefore, to obtain a predicted distribution of the same functional form as the posterior \eqref{eq:multiple_ellipse_ET_distribution}, following the discussion in \cite{GranstromO:2014} we solve independent integrals instead,
\begin{subequations}
\begin{align}
p\left(\sx_{+}\left|\setZ^{k}\right.\right) = & \int p\left(\sx_{+} \left| \sx_{} \right.\right) p\left(\sx_{}\left|\setZ^{k}\right.\right)  \diff\sx_{}, \label{eq:kinematics_prediction_integral} \\
p\left(\ext_{+}^{(i)} \left| \setZ^{k} \right.\right) = & \int p\left(\ext_{+}^{(i)} \left| \sx_{} , \setZ^{k} \right. \right) p\left(\sx_{}\left|\setZ^{k}\right.\right)  \diff\sx_{}. \label{eq:extension_prediction_integral}
\end{align}%
\end{subequations}%
%\begin{align}
%& \left( \prod_{i=1}^{N_{s,k}} \int p\left(\ext_{+}^{(i)} \left| \sx_{} , \setZ^{k} \right.\right) p\left(\sx_{}\left|\setZ^{k}\right.\right)  \diff\sx_{} \right) \nonumber \\
%& \times \int p\left(\sx_{+} \left| \sx_{} \right.\right) p\left(\sx_{}\left|\setZ^{k}\right.\right)  \diff\sx_{} \\
%= & p\left(\sx_{+}\left|\setZ^{k}\right.\right) \left( \prod_{i=1}^{N_{s,k}} p\left(\ext_{+}^{(i)} \left| \setZ^{k} \right.\right) \right). \nonumber
%\end{align}
With a mixture distribution, the predicted mixture is the mixture of the predicted components,
\begin{subequations}
\begin{align}
p\left(\xi_{+}\left|\setZ^{k}\right.\right) = & \int p(\xi_{+}|\xi) \sum_{\ell=1}^{J} w^{(\ell)} p^{(\ell)} \left(\xi \left| \setZ^{k}\right.\right) \diff \xi_{} \\
%= & \sum_{\ell=1}^{J} w^{(\ell)} \int p(\xi_{+}|\xi_k) p^{(\ell)} \left(\xi_{k} \left| \setZ^{k}\right.\right) \diff \xi_{} \\
= & \sum_{\ell=1}^{J} w^{(\ell)} p^{(\ell)}\left(\xi_{+}\left|\setZ^{k}\right.\right).
\end{align}
\end{subequations}

\subsection{Correction}
\label{sec:multi_ellipse_correction}
Let $\theta$ denote a possible measurement-to-subobject association event, and let $\setTheta$ denote the set of all possible association events.  For measurement generation, we assume the following:
\begin{assumption}
The subobjects generate measurements independently of each other.
For each subobject, the generated measurements are independent.
Each measurement is generated by exactly one subobject.
Measurement origin is unknown. \hfill$\square$
\end{assumption}
\begin{remark}
These assumptions are analogous to multiple target tracking, where it is typically assumed that each target generates measurements independently of the other targets, that the target generated measurements are independent, that each measurement is generated by exactly one target, and that measurement origin is unknown, see \egp \cite{BarShalomWT:2011}.
\hfill$\square$\end{remark}

Under an association event $\theta$ the measurement set $\setZ_{k}$ can be partitioned into $N_{s,k}$ (possibly empty) subsets,
\begin{align}
\setZ_{k} = & \bigcup_{i=1}^{N_{s,k}} \setZ_{k}^{(\theta,i)}, & \setZ_{k}^{(\theta,i)} = & \left\{\sz_{k}^{(\theta,i,j)}\right\}_{j=1}^{n_{z,k}^{(\theta,i)}},
\end{align}
where the $i$th subset $\setZ_{k}^{(\theta,i)}$ was generated by the $i$th subobject. Conditioned on $\theta$ the measurement likelihood is
\begin{subequations}
\begin{align}
p\left(\setZ_{k}\left|\xi_{k},\theta\right.\right) = & \prod_{i=1}^{N_{s,k}} p\left(\setZ_{k}^{(\theta,i)}\left|\gamma_{k}^{(i)},\sx_{k},\ext_{k}^{(i)}\right.\right).
\end{align}%
\end{subequations}%
If the $i$th subset is empty (\iep $n_{z,k}^{(\theta,i)}=0$) the subset likelihood is simply the likelihood of an empty set of measurements,
\begin{align}
& p\left(\setZ_{k}^{(\theta,i)}\left|\gamma_{k}^{(i)},\sx_{k},\ext_{k}^{(i)}\right.\right) = P\left(n_{z,k}^{(\theta,i)}=0\left| \gamma_{k}^{(i)}\right.\right).
\end{align}%
If $n_{z,k}^{(\theta,i)}>0$ the subobject likelihood is
\begin{subequations}
\begin{align}
p&  \left(\setZ_{k}^{(\theta,i)}\left|\gamma_{k}^{(i)},\sx_{k},\ext_{k}^{(i)}\right.\right) \nonumber \\
%= & n_{z,k}^{(\theta,i)}! P\left(n_{z,k}^{(\theta,i)}\left| \gamma_{k}^{(i)},\sx_{k},\ext_{k}^{(i)} \right.\right) \nonumber \\
%& \times p\left(\setZ_{k}^{(\theta,i)}\left|\gamma_{k}^{(i)},\sx_{k},\ext_{k}^{(i)},n_{z,k}^{(\theta,i)}\right.\right) \\
%= & n_{z,k}^{(\theta,i)}! P\left(n_{z,k}^{(\theta,i)}\left| \gamma_{k}^{(i)}\right.\right) p\left(\setZ_{k}^{(\theta,i)}\left|\sx_{k},\ext_{k}^{(i)}\right.\right) \\
= & n_{z,k}^{(\theta,i)}! P\left(n_{z,k}^{(\theta,i)}\left| \gamma_{k}^{(i)}\right.\right) \prod_{j=1}^{n_{z,k}^{(\theta,i)}} p\left(\sz_{k}^{(\theta,i,j)}\left|\sx_{k},\ext_{k}^{(i)}\right.\right).
\end{align}%
\end{subequations}%
Let the predicted mixture distribution be
\begin{align}
p \left(\xi_{k} \left| \setZ^{k-1}\right.\right) = \sum_{\ell=1}^{J} w^{(\ell)} p^{(\ell)} \left(\xi_{k} \left| \setZ^{k-1}\right.\right).
\end{align}
By the total probability theorem the density $p\left(\xi_{k}|\setZ^{k}\right)$ is %\cite{LanRL:2014}
\begin{align}
p\left(\xi_{k}|\setZ^{k}\right) = & \sum_{\theta\in\setTheta} p\left(\xi_{k}|\setZ^{k},\theta\right)P\left(\theta|\setZ^{k}\right),
\label{eq:Bayes_updated_all_assoc_events}
\end{align}
where $p\left(\xi_{k}|\setZ^{k},\theta\right)$ is the Bayes updated distribution for the association event $\theta$, and $P\left(\theta|\setZ^{k}\right)$ is the probability of the association event $\theta$. As noted in \cite[Assumption 3]{LanRL:2014}, without any prior information the association events can be assumed to be equally likely, \iep $P\left(\theta|\setZ^{k-1}\right)= |\setTheta|^{-1}$. In this case we have
\begin{subequations}
\begin{align}
P\left(\theta|\setZ^{k}\right) = & \frac{p\left(\setZ_{k}|\theta,\setZ^{k-1}\right)P\left(\theta|\setZ^{k-1}\right)}{\sum_{\theta'\in\setTheta}p\left(\setZ_{k}|\theta',\setZ^{k-1}\right)P\left(\theta'|\setZ^{k-1}\right)} \\
%= & \frac{p\left(\setZ_{k}|\theta,\setZ^{k-1}\right)}{\sum_{\theta'\in\setTheta}p\left(\setZ_{k}|\theta',\setZ^{k-1}\right)} \\
= & \frac{\sum_{\ell=1}^{J} w^{(\ell)} p^{(\ell)}\left(\setZ_{k}|\theta,\setZ^{k-1}\right)}{\sum_{\theta'\in\setTheta} \sum_{\ell'=1}^{J} w^{(\ell')} p^{(\ell')}\left(\setZ_{k}|\theta',\setZ^{k-1}\right)},
\end{align}%
\label{eq:prob_assoc_event}%
\end{subequations}%
where we have again used the total probability theorem in the second equality. For the association event $\theta$ the Bayes updated distribution is
\begin{subequations}
\begin{align}
p\left(\xi_{k}\left|\setZ^{k},\theta\right.\right) = & \frac{p\left(\setZ_{k}\left|\xi_{k},\theta\right.\right) p\left(\xi_{k}\left|\setZ^{k-1}\right.\right)}{\int p\left(\setZ_{k}\left|\xi_{k},\theta\right.\right) p\left(\xi_{k}\left|\setZ^{k-1}\right.\right) \diff \xi_{k}} \\
%= & \frac{p\left(\setZ_{k}\left|\xi_{k},\theta\right.\right) \sum_{\ell=1}^{J} w^{(\ell)} p^{(\ell)} \left(\xi_{k} \left| \setZ^{k-1}\right.\right)}{\int p\left(\setZ_{k}\left|\xi_{k},\theta\right.\right) \sum_{\ell=1}^{J} w^{(\ell)} p^{(\ell)} \left(\xi_{k} \left| \setZ^{k-1}\right.\right) \diff \xi_{k}} \\
%= & \frac{ \sum_{\ell=1}^{J} w^{(\ell)} p\left(\setZ_{k}\left|\xi_{k},\theta\right.\right) p^{(\ell)} \left(\xi_{k} \left| \setZ^{k-1}\right.\right)}{\sum_{\ell=1}^{J} w^{(\ell)} \int p\left(\setZ_{k}\left|\xi_{k},\theta\right.\right) p^{(\ell)} \left(\xi_{k} \left| \setZ^{k-1}\right.\right) \diff \xi_{k}} \\
= & \frac{ \sum_{\ell=1}^{J} w^{(\ell)} p^{(\ell)}\left(\setZ_{k}\left|\theta,\setZ^{k-1}\right.\right) p^{(\ell)} \left(\xi_{k} \left| \setZ^{k},\theta\right.\right)}{\sum_{\ell=1}^{J} w^{(\ell)} p^{(\ell)}\left(\setZ_{k}\left|\theta,\setZ^{k-1}\right.\right) }.
%= & \frac{p\left(\setZ_{k}\left|\xi_{k},\theta\right.\right) p\left(\xi_{k}\left|\setZ^{k-1}\right.\right)}{p\left(\setZ_{k}\left|\setZ^{k-1},\theta\right.\right)}
\end{align}%
\label{eq:Bayes_updated_single_assoc_event}%
\end{subequations}%
Combining \eqref{eq:Bayes_updated_all_assoc_events}, \eqref{eq:prob_assoc_event} and \eqref{eq:Bayes_updated_single_assoc_event} gives the posterior distribution
\begin{align}
p\left(\xi_{k}|\setZ^{k}\right) = & \sum_{\theta\in\setTheta} \sum_{\ell=1}^{J} w^{(\ell)}\left(\theta\right) p^{(\ell)} \left(\xi_{k} \left| \setZ^{k},\theta\right.\right), \\
w^{(\ell)}\left(\theta\right) = & \frac{w^{(\ell)} p^{(\ell)}\left(\setZ_{k}|\theta,\setZ^{k-1}\right)}{\sum_{\theta'\in\setTheta} \sum_{\ell'=1}^{J} w^{(\ell')} p^{(\ell')}\left(\setZ_{k}|\theta',\setZ^{k-1}\right)},
\end{align}
where, following the assumption that the subobjects generate measurements independently, for the predicted likelihood we have
\begin{subequations}
\begin{align}
& p^{(\ell)}\left(\setZ_{k}|\theta,\setZ^{k-1}\right) = \prod_{i=1}^{N_{s,k}} p^{(\ell)}\left(\left.\setZ_{k}^{(\theta,i)}\right|\setZ^{k-1}\right). \label{eq:meas_lik_prod} \\
& p\left(\setZ_{k}|\setZ^{k-1}\right) = \frac{1}{|\setTheta|} \sum_{\theta\in\setTheta} \sum_{\ell=1}^{J} w^{(\ell)} p^{(\ell)}\left(\setZ_{k}|\theta,\setZ^{k-1}\right),
\end{align}
\end{subequations}
The predicted likelihood $p\left(\setZ_{k}|\setZ^{k-1}\right)$ is useful in a multiple target tracking scenario, \egp if the presented extended target model is used in an implementation of an extended target \phd or \cphd filter \cite{mahler_FUSION_2009_extTarg,LundquistGO:2013}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%
%%%%% GGIW implementation
%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{A Gamma Gaussian inverse Wishart implementation}
\label{sec:GGIW_implementation}
In this section we give a gamma Gaussian inverse Wishart implementation of the multiple random matrix extended target model outlined above. To handle different types of motion $M_{k}$ different motion models are used.

\subsection{Extended target state distribution}
In Gilholm \etal's extended target model \cite{GilholmGMS:2005,GilholmS:2005} the number of measurements that each target generates is Poisson distributed with a parameter $\varsigma\left(\xi\right)$ that is a function of the extended target state. In practice this means that the expected value of the number of measurements generated by a target with state $\xi$ is $\varsigma\left(\xi\right)$. Here the Poisson model is adopted for each subobject and the parameters are given by the measurement rates, \iep $\varsigma^{(i)}\left(\xi_{k}\right) = \gamma_{k}^{(i)}, \ \forall i$.
%\begin{align}
%\varsigma^{(i)}\left(\xi_{k}\right) = \gamma_{k}^{(i)}, \quad \forall i.
%\end{align}
The gamma distribution is the conjugate prior for the Poisson distribution's parameter, and the subobjects' measurement rates are modeled as gamma distributed.
%\cite{GranstromO:2012c}
%\begin{align}
%p\left(\gamma_{k}^{(i)} \left| \setZ^{k}\right.\right) = \Gammapdf{\gamma_{k}^{(i)}}{\alpha_{k|k}^{(i)}}{\beta_{k|k}^{(i)}}.
%\end{align}
Following the random matrix model \cite{Koch:2008,FeldmannFK:2011} the subobjects' random matrices are modeled as inverse Wishart distributed.
%\begin{align}
%p\left(\ext_{k}^{(i)} \left| \setZ^{k} \right.\right) = \IWishpdf{\ext_{k}^{(i)}}{v_{k|k}^{(i)}}{V_{k|k}^{(i)}}.
%\end{align}
The kinematic vector is modeled as Gaussian distributed.
%\begin{align}
%p\left(\sx_{k}\left|\setZ^{k}\right.\right) = \Npdfbig{\sx_{k}}{m_{k|k}}{P_{k|k}}.
%\end{align}
The extended target state distribution is
\begin{subequations}
\begin{align}
p\left(\xi_{k} \left| \setZ^{k}\right.\right) = & \Npdfbig{\sx_{k}}{m_{k|k}}{P_{k|k}} \prod_{i=1}^{N_{s,k}} \bigg( \Gammapdf{\gamma_{k}^{(i)}}{\alpha_{k|k}^{(i)}}{\beta_{k|k}^{(i)}} \big. \nonumber \\
& \times \big. \IWishpdf{\ext_{k}^{(i)}}{v_{k|k}^{(i)}}{V_{k|k}^{(i)}} \bigg) \label{eq:GGIW_extended_target_state} \\
= & \GGIWpdf{\xi_{k}}{\zeta_{k|k}},
\end{align}
\end{subequations}
where $\GGIWpdf{\cdot}{\cdot}$ is introduced for brevity and $\zeta_{k|k}$ denotes all of the involved parameters. As noted above, the extended target state distribution is represented by a distribution mixture
\begin{align}
p\left(\xi_{k} \left| \setZ^{k}\right.\right) = & \sum_{\ell=1}^{J_{k|k}} w_{k|k}^{(\ell)} \GGIWpdf{\xi_{k}}{\zeta_{k|k}^{(\ell)}}, \label{eq:mixtureGGIW_extended_target_state}
\end{align}
where $\sum_{\ell}w_{k|k}^{(\ell)}=1$.

%\subsection{Target estimate initialization}
When a new target appears the parameters $\zeta^{(\ell)}$ of the estimate must be initialized. Table~\ref{tab:multi_ellipse_initilization} gives a simple algorithm where this is performed using the first set of measurements. The algorithm initializes $N_{p}$ hypotheses in each motion mode.
A simple initialization example is given in Figure~\ref{fig:Initialization_Example_1}. In this example there is a single motion mode, and $N_{p}=4$ hypotheses are generated using only $8$ measurements.

\begin{figure}[!t]
\centerline{
\subfloat{\includegraphics[width=0.25\columnwidth]{Initialization_Example_1} \label{fig:Initialization_Example_1}}
\subfloat{\includegraphics[width=0.25\columnwidth]{Initialization_Example_2} \label{fig:Initialization_Example_2}}
\subfloat{\includegraphics[width=0.25\columnwidth]{Initialization_Example_3} \label{fig:Initialization_Example_3}}
\subfloat{\includegraphics[width=0.25\columnwidth]{Initialization_Example_4} \label{fig:Initialization_Example_4}}
}
\caption{Initialization example. True underlying extended target (orange area), measurements (red squares), and initialized estimates (blue ellipses, initialized main subobject shown by thicker line).}
\label{fig:Initialization_Example}
\end{figure}


\begin{table}[tb]
\caption{Multiple ellipse parameter initialization}
\label{tab:multi_ellipse_initilization}
%\centering
\vspace{-0.5cm}
\rule[0pt]{\columnwidth}{1pt}
\vspace{-0.3cm}
\begin{algorithmic}[1]
\STATE \textbf{Input:} Set of measurements $\setZ=\{\sz_{i}\}_{i=1}^{n}$. Desired number of initial hypotheses $N_{p}$. Initial kinematics $\mathbf{c}_{0}$ and initial covariance $P_{0}$. Initial mean $e$ and variance $v$ for measurement rates.
\STATE Define $\sz_{c} = \frac{1}{n}\sum_{i=1}^{n}\sz_{i} $, $r_{z} = \frac{1}{2}\max_{i} \left\| \sz_{i}-\sz_{c} \right\|_{2}$. Set $\ell=0$.
\FOR{$p=1,\ldots,N_{p}$}
	\FOR{$m=1,\ldots,M_{k}$}
		\STATE Set $\ell = \ell+1$
		\STATE $\gamma$: $\alpha_{0}^{(\ell,i)} = \frac{e^2}{v}$, $\beta_{0}^{(\ell,i)} = \frac{e}{v}$.
		\STATE $\sx$: $P_{0}^{(\ell)}=P_{0}$, $\mathbf{p}_{0}^{(\ell)}=\sz_{c}$, $\mathbf{c}_{0}^{(\ell)}=\mathbf{c}_{0}$,
			\begin{align*}
				\mathbf{d}_{0}^{(\ell,i)} = r_{z}\begin{bmatrix} \cos\left(\frac{2\pi(i-2)}{N_{s}-1}+\frac{2\pi(p-1)}{N_{s}N_{p}}\right) \\ \sin\left(\frac{2\pi(i-2)}{N_{s}-1}+\frac{2\pi(p-1)}{N_{s}N_{p}}\right) \end{bmatrix}.
			\end{align*}
		\STATE $\ext$: $v_{0}^{(\ell,i)} = 2d+5$, $V_{0}^{(\ell,i)} = \left(\frac{r_z}{4}\right)^{2}\Id\left(v_{0}^{(\ell,i)}-2d-2\right)$.
	\ENDFOR
\ENDFOR
\STATE \textbf{Output:} $p\left(\xi_{0}\right) = \sum_{\ell=1}^{J_{0}} w_{0}^{(\ell)} \GGIWpdf{\xi_{0}}{\zeta_{0}^{(\ell)}}$ where $w_{0}^{(\ell)} = \frac{1}{J_{0}}$.
\end{algorithmic}
\vspace{-0.1cm}
\rule[0pt]{\columnwidth}{1pt}
\end{table}

\subsection{Prediction}
With a posterior distribution of the form \eqref{eq:mixtureGGIW_extended_target_state} the predicted distribution is 
%\begin{subequations}
\begin{align}
%p\left(\xi_{k+1} \left| \setZ^{k}\right.\right) = & \left( \prod_{i=1}^{N_{s,k}} \Gammapdf{\gamma_{k+1}^{(i)}}{\alpha_{k+1|k}^{(i)}}{\beta_{k+1|k}^{(i)}} \right) \nonumber \\
%& \times \Npdfbig{\sx_{k+1}}{m_{k+1|k}^{(\theta)}}{P_{k+1|k}^{(\theta)}} \nonumber \\
%& \times \left( \prod_{i=1}^{N_{s,k}} \IWishpdf{\ext_{k+1}^{(i)}}{v_{k+1|k}^{(i)}}{V_{k+1|k}^{(i)}} \right) 
p\left(\xi_{k+1} \left| \setZ^{k}\right.\right) = & \sum_{m=1}^{M_{k}}\sum_{\ell=1}^{J_{k|k}} \pi_{m,m'(\ell)} w_{k|k}^{(\ell)} \GGIWpdf{\xi_{k}}{\zeta_{k+1|k}^{(m,\ell)}}, \label{eq:predicted_GGIW_distribution}%
%= & \sum_{\ell=1}^{J_{k+1|k}}w_{k+1|k}^{(\ell)}\GGIWpdf{\xi_{k}}{\zeta_{k+1|k}^{(\ell)}},
\end{align}%
%\end{subequations}%
where $\pi_{m,m'(\ell)}$ is the probability of a transition to the current mode $m$ from the previous mode $m'(\ell)$ that component $\ell$ was in.

\subsubsection{Measurement rates}
For the measurement rates the exponential forgetting prediction from \cite{GranstromO:2012c} is used. For the $m$th motion model the parameters are predicted as
%\begin{subequations}
\begin{align}
\alpha_{k+1|k}^{(m,\ell,i)} & = \frac{\alpha_{k|k}^{(\ell,i)}}{\eta_k^{(m)}}, \qquad \beta_{k+1|k}^{(m,\ell,i)} =  \frac{\beta_{k|k}^{(\ell,i)}}{\eta_k^{(m)}},
\end{align}%
%\end{subequations}%
which corresponds to keeping the expected value of $\gamma_{k}^{(i)}$ constant, while increasing the variance with a factor $\eta_{k}^{(m)}$ \cite{GranstromO:2012c}. This prediction has an effective window length of $w_{e} =\frac{\eta_{k}^{(m)}}{\eta_{k}^{(m)}-1}$, where $\frac{1}{\eta_k^{(m)}} < 1$ is the forgetting factor.

\subsubsection{Kinematic state}
For the $m$th motion model the kinematic state transition density is modeled as
\begin{align}
p(\sx_{k+1}|\sx_k) = & \Npdfbig{\sx_{k+1}}{f^{(m)}(\sx_k)}{Q_{k+1}^{(m)}},
\end{align}
where $f^{(m)}(\cdot):\mathbb{R}^{n_x}\rightarrow\mathbb{R}^{n_x}$ is a state transition function, and $Q_{k+1}^{(m)}$ is the process noise covariance for the kinematic state. The transition function can be partitioned into $N_{s,k}$ parts,
\begin{align}
%f^{(m)}(\sx_k) = & \begin{bmatrix} f_{\mathbf{p},\mathbf{c}}^{(m)}\left(\mathbf{p}_{k},\mathbf{c}_{k}\right) \\ f_{\mathbf{d}}^{(m)}\left(\mathbf{d}_{k}^{(2)},\mathbf{c}_{k}\right) \\ \vdots \\ f_{\mathbf{d}}^{(m)}\left(\mathbf{d}_{k}^{(N_{s,k})},\mathbf{c}_{k}\right)\end{bmatrix}, \\
f^{(m)}(\sx_k) = & \begin{bmatrix} f_{\mathbf{p},\mathbf{c}}^{(m)}\left(\mathbf{p}_{k},\mathbf{c}_{k}\right)^{\tp} & \ldots & f_{\mathbf{d}}^{(m)}\left(\mathbf{d}_{k}^{(i)},\mathbf{c}_{k}\right)^{\tp} & \ldots \end{bmatrix}^{\tp},
\end{align}
%\begin{align}
%f^{(m)}(\sx_k) = & \begin{bmatrix} f_{\mathbf{p},\mathbf{c}}^{(m)}\left(\mathbf{p}_{k},\mathbf{c}_{k}\right) \\ f_{\mathbf{d}}^{(m)}\left(\mathbf{d}_{k}^{(2)},\mathbf{c}_{k}\right) \\ \vdots \\ f_{\mathbf{d}}^{(m)}\left(\mathbf{d}_{k}^{(i)},\mathbf{c}_{k}\right) \\ \vdots \\ f_{\mathbf{d}}^{(m)}\left(\mathbf{d}_{k}^{(N_{s,k})},\mathbf{c}_{k}\right)\end{bmatrix},
%\end{align}
where $f_{\mathbf{p},\mathbf{c}}^{(m)}(\cdot)$ describes the time evolution of the overall position and the kinematics, and $f_{\mathbf{d}}^{(m)}(\cdot)$ describes the time evolution of the subobject offsets. The functions $f_{\mathbf{p},\mathbf{c}}^{(m)}(\cdot)$ and $f_{\mathbf{d}}^{(m)}(\cdot)$ are generally nonlinear, see \cite{RongLiJ:2003} for a thorough overview of state transition functions.

In case $f^{(m)}\left(\sx_{k}\right)$ is a linear function, the solution to \eqref{eq:kinematics_prediction_integral}
%\begin{align}
%p\left(\sx_{k+1}\left|\setZ^{k}\right.\right) =\int & \Npdfbig{\sx_{k+1}}{f^{(m)}(\sx_k)}{Q_{k+1}^{(m)}} \nonumber \\
%& \times \Npdfbig{\sx_k}{\mkk}{P_{k|k}}\diff\sx_k. \label{eqn:predictionkinematic}
%\end{align}
is given by the Kalman filter prediction \cite{Kalman:1960}. If $f^{(m)}(\sx_{k})$ is non-linear it is straightforward to solve \eqref{eq:kinematics_prediction_integral} approximately. Using the extended Kalman filter prediction formulas, see \egp \cite{Jazwinski:1970}, the predicted mean $m_{k+1|k}^{(m,\ell)}$ and covariance $P_{k+1|k}^{(m,\ell)}$ are
\begin{subequations}
\begin{align}
m_{k+1|k}^{(m,\ell)} & = f^{(m)}(\mkk^{(\ell)}),\\
P_{k+1|k}^{(m,\ell)} & = F_{k|k}^{(m,\ell)}P_{k|k}^{(\ell)}\left(F_{k|k}^{(m,\ell)}\right)^{\tp}+Q_{k+1}^{(m)}
\end{align}%
\end{subequations}%
where $F_{k|k}^{(m,\ell)}=\left.\nabla_{\sx} f^{(m)}(\sx)\right|_{\sx=\mkk^{(\ell)}}$ is the gradient of $f^{(m)}(\cdot)$ evaluated at the mean $\mkk^{(\ell)}$.

\subsubsection{Random matrices}
For the $m$th motion model we use the transition density suggested in \cite{GranstromO:2014},
\begin{align}
p&(\ext_{k+1}^{(i)}|\sx_{k},\ext_{k}^{(i)}) \\
% & = \Wishpdf{\ext_{k+1}^{(i)}}{n_{k+1}^{(m)}}{\frac{M_{\sx_{k}}^{(m)}\ext_{k}^{(i)}\left(M^{(m)}_{\sx_{k}}\right)^{\tp}}{n_{k+1}^{(m)}}},
& = \Wishpdf{\ext_{k+1}^{(i)}}{n_{k+1}^{(m)}}{\left(n_{k+1}^{(m)}\right)^{-1}M_{\sx_{k}}^{(m)}\ext_{k}^{(i)}\left(M^{(m)}_{\sx_{k}}\right)^{\tp}}, \nonumber
\end{align}
where $n_{k+1}^{(m)}>d-1$ is a scalar design parameter and the matrix transformation $M_{\sx_k}^{(m)}\triangleq M^{(m)}\left(\sx_k\right):\mathbb{R}^{n_x}\rightarrow \mathbb{R}^{d\times d}$ is a non-singular matrix valued function of the kinematic state. The extension state's time evolution is modeled as being dependent on the kinematic state mainly because it allows for the modeling of rotation of extended targets, however in general the only requirement is that the output is a non-singular $d\times d$ matrix \cite{GranstromO:2014}.
%For the random matrices we have
%\begin{align}
%& \iint p\left(\ext_{k+1}^{(i)} \left| \sx_{k},\ext_{k}^{(i)} \right.\right) p\left(\sx_{k}\left|\setZ^{k}\right.\right) p\left(\ext_{k}^{(i)} \left| \setZ^{k}\right.\right) \diff\ext_{k}^{(i)} \diff\sx_{k}.
%\end{align}
Details on how the parameters $v_{k+1|k}^{(m,\ell,i)}$ and $V_{k+1|k}^{(m,\ell,i)}$ of the solutions to \eqref{eq:extension_prediction_integral}
%predicted distributions
%\begin{align}
%p\left(\ext_{k+1}^{(i)}|\setZ^{k}\right)= & \IWishpdf{\ext_{k+1}^{(i)}}{v_{k+1|k}^{(m,\ell,i)}}{V_{k+1|k}^{(m,\ell,i)}}
%\end{align}
are computed are given in \cite{GranstromO:2014}.

\subsection{Generation of association events}
\label{sec:subset_association_events}
The correction step, see Section~\ref{sec:multi_ellipse_correction}, involves a summation over the set $\setTheta$ of all possible measurement-to-subobject association events. For $n_{z,k}$ measurements and $N_{s,k}$ subobjects there are $(N_{s,k})^{n_{z,k}}$ possible measurement-to-subobject association events. For example, if $N_{s,k}=3$ and $n_{z,k}=5$ there are $243$ possible association events, and if $n_{z,k}=10$ there are $59049$ possible association events. Due to the quickly increasing size of the full set of association events approximations are necessary to achieve tractable computational complexity.

Using different methods to simplify the data association problem is common in target tracking. Popular methods for multiple point target tracking include probabilistic data association (\textsc{pda}), and multiple hypothesis tracking (\textsc{mht}), see e.g. \cite{BarShalomWT:2011}. Data clustering methods are used in extended target \phd/\cphd filters to reduce the number of measurement partitions that are considered, see e.g. \cite{GranstromLO:2010,GranstromLO:2012,GranstromO:2012a,LundquistGO:2013}. In \cite{LanRL:2014} it is proposed to reduce the number of measurement-to-subobject association events by using a combination of $k$-means clustering, see e.g. \cite{Bishop:2006,HastieTF:2009}, and gating with a suitable pseudo-likelihood.

%This is similar to the extended target \phd and \cphd filters \cite{mahler_FUSION_2009_extTarg,LundquistGO:2013}, which both require a summation over all possible partitions of the measurement set. The number of possible partitions of a set increases rapidly when the size of the set increases, and in the implementations \cite{GranstromLO:2010,GranstromLO:2012,GranstromO:2012a,LundquistGO:2013} the full set of partitions is approximated with a subset to achieve computational tractability.

In this paper a subset $\bar{\setTheta}\subseteq\setTheta$ of association events is computed using a method that is based on the Expectation Maximization algorithm \cite{DempsterLR:1977} for Gaussian Mixtures (\emgm), see \egp \cite[Chapter 9]{Bishop:2006}. \emgm is favored over $k$-means clustering because it has been shown that the $k$-means clustering algorithm often can give unfavorable results for elliptically shaped extended targets, see \cite{GranstromO:2012a}. %, see Table~\ref{tab:generation_of_association_events}.
First \emgm is used to partition the current set of measurements into $N_{c}$ clusters. To accommodate the possibility that one or more subobjects does not generate any measurements \emgm is used for $N_{c}\in[1, \ 2,\ldots,N_{s,k}]$. Because the solution to \emgm has many local optimums, for each $N_{c}$ the algorithm is given several different initializations. Note that care is taken to ensure that the set of partitions returned by \emgm only contains unique partitions.

The next step is to use the clusters to obtain measurement-to-subobject associations. Given a partition of the measurement set with $N_{c}$ clusters, and an estimate with $N_{s,k}$ subobjects, there are ${N_{s,k}!}/{(N_{s,k}-N_{c})!}$
%\begin{align}
%\frac{N_{s,k}!}{(N_{s,k}-N_{c})!}
%\end{align}
possible cluster-to-subobject associations. A cluster-to-subobject association defines a measurement-to-subobject association event $\theta$ because each measurement is associated to a cluster, which in turn is associated to a subobject. Let $C(N_{c})$ denote the number of unique partitions with $N_{c}$ clusters obtained using \emgm. Then the number of measurement-to-subobject association events that has to be considered is
\begin{align}
\left|\bar{\setTheta}\right| = \sum_{N_c=1}^{N_{s,k}} C(N_{c}) \frac{N_{s,k}!}{\left(N_{s,k}-N_{c}\right)!}.
\end{align}
% Empirically we have found that this number typically is several orders of magnitude smaller than $\left(N_{s,k}\right)^{n_{z,k}}$.

\subsection{Correction}
Each object generates a Poisson distributed number of measurements and the measurement models are linear Gaussian,
\begin{subequations}
\begin{align}
P\left(n_{z,k}^{(i)}\left| \gamma_{k}^{(i)}\right.\right) = & \mathcal{PS}\left(n_{z,k}^{(i)};\ \gamma_{k}^{(i)}\right) \\
p\left(\sz_{k}^{(i,j)}\left|\sx_{k},\ext_{k}^{(i)}\right.\right) = & \Npdfbig{\sz_{k}^{(i,j)}}{H_{k}^{(i)}\sx_{k}}{\ext_{k}^{(i)}}. \label{eq:LinGaussMeasurementModel}
\end{align}
\end{subequations}
For the kinematic state \eqref{eq:kinematic_state_vector} the models $H_{k}^{(i)}$ are
\begin{subequations}
\begin{align}
H_{k}^{(1)} = & \begin{bmatrix} \Id & \mathbf{0}_{d\times n_c} & \mathbf{0}_{d\times (N_{s,k}-1)d} \end{bmatrix},\\
H_{k}^{(i)} = & \begin{bmatrix} \Id & \mathbf{0}_{d\times n_c} & \mathbf{0}_{d\times (i-2)d} & \Id & \mathbf{0}_{d\times (N_{s,k}-i)d} \end{bmatrix},
\end{align}%
\end{subequations}%
for $i=2,\ldots,N_{s,k}$.
% Let $\bar{\setTheta} \subseteq \setTheta$ be a subset of all possible association events.
For an association event $\theta\in\bar{\setTheta}$ the centroid measurement and scatter matrix are defined as follows,
\begin{subequations}
\begin{align}
\bar{\sz}_{k}^{(\theta,i)} & = \frac{1}{n_{z,k}^{(\theta,i)}} \sum_{j=1}^{n_{z,k}^{(\theta,i)}} \sz_{k}^{(\theta,i,j)},\\
Z_{k}^{(\theta,i)} & = \sum_{j=1}^{n_{z,k}^{(\theta,i)}} \left(\sz_{k}^{(\theta,i,j)}-\bar{\sz}_{k}^{(\theta,i)}\right) \left(\sz_{k}^{(\theta,i,j)}-\bar{\sz}_{k}^{(\theta,i)}\right)^{\tp}.
\end{align}
\end{subequations}
The same measurement model is used for all motion models. With a predicted distribution %of the form \eqref{eq:predicted_GGIW_distribution}
\begin{align}
p\left(\xi_{k} \left| \setZ^{k-1}\right.\right) = & \sum_{\ell=1}^{J_{k|k-1}} w_{k|k-1}^{(\ell)} \GGIWpdf{\xi_{k}}{\zeta_{k|k-1}^{(\ell)}}
\end{align}
the corrected distribution is 
\begin{align}
p\left(\xi_{k} \left| \setZ^{k}\right.\right) = & \sum_{\theta\in\bar{\setTheta}}\sum_{\ell=1}^{J_{k|k-1}}w_{k|k}^{(\theta,\ell)}\GGIWpdf{\xi_{k}}{\zeta_{k|k}^{(\theta,\ell)}}
%= & \sum_{\ell=1}^{J_{k|k}}w_{k|k}^{(\ell)}\GGIWpdf{\xi_{k}}{\zeta_{k|k}^{(\ell)}} \label{eq:corrected_GGIW_distribution}
%\sum_{\theta\in\bar{\setTheta}} w_{k|k}^{\theta} \left( \prod_{i=1}^{N_{s,k}} \Gammapdf{\gamma_{k}^{(i)}}{\alpha_{k|k}^{(\theta,i)}}{\beta_{k|k}^{(\theta,i)}} \right) \nonumber \\
%& \times \Npdfbig{\sx_{k}}{m_{k|k}^{(\theta)}}{P_{k|k}^{(\theta)}} \nonumber \\
%& \times \left( \prod_{i=1}^{N_{s,k}} \IWishpdf{\ext_{k}^{(i)}}{v_{k|k}^{(\theta,i)}}{V_{k|k}^{(\theta,i)}} \right). \label{eq:corrected_GGIW_distribution}
\end{align}
In the sections that follow the parameters $\zeta_{k|k}^{(\theta,\ell)}$ of the corrected distribution are given. The derivation of the measurement update is given in the Appendix. %\cite{Granstrom:2014proof}.

% A derivation of the measurement update is given in the Appendix.%~\ref{app:measurement_update}.

\subsubsection{Measurement rates}
The corrected parameters are
%\begin{subequations}
\begin{align}
\alpha_{k|k}^{(\theta,\ell,i)} & = \alpha_{k|k-1}^{(\ell,i)} + n_{z,k}^{(\theta,i)}, \qquad \beta_{k|k}^{(\theta,\ell,i)} = \beta_{k|k-1}^{(\ell,i)} + 1.
\end{align}
%\end{subequations}
%where $n_{z,k}^{(\theta,i)}$ are the number of measurements associated to subobject $i$ under the association event $\theta$.

\subsubsection{Kinematic state}
The corrected parameters are
\begin{subequations}
\begin{align}
m_{k|k}^{(\theta,\ell)} & = m_{k|k-1}^{(\ell)} + K_{k}^{(\theta,\ell)}\left(\bar{\sz}_{k}^{(\theta)}-\mathbb{H}_{k}m_{k|k-1}^{(\ell)}\right),\\
P_{k|k}^{(\theta,\ell)} & = P_{k|k-1}^{(\ell)} + K_{k}^{(\theta,\ell)}\mathbb{H}_{k}P_{k|k-1}^{(\ell)},\\
\bar{\sz}_{k}^{(\theta)} & = \begin{bmatrix} \left(\bar{\sz}_{k}^{(\theta,1)}\right)^{\tp} & \cdots & \left(\bar{\sz}_{k}^{(\theta,N_{s,k})}\right)^{\tp}  \end{bmatrix}^{\tp}, \\
\mathbb{H}_{k} & = \begin{bmatrix} \left(H_{k}^{(1)}\right)^{\tp} & \cdots & \left(H_{k}^{(N_{s,k})}\right)^{\tp} \end{bmatrix}^{\tp},\\
K_{k}^{(\theta,\ell)} & = P_{k|k-1}^{(\ell)}\mathbb{H}_{k}^{\tp}\left(S_{k}^{(\theta,\ell)}\right)^{-1},\\
S_{k}^{(\theta,\ell)} & = \mathbb{H}_{k}P_{k|k-1}^{(\ell)}\mathbb{H}_{k}^{\tp} + \hat{\mathbb{X}}_{k|k-1}^{(\theta,\ell)},\\
\hat{\mathbb{X}}_{k|k-1}^{(\theta,\ell)} & = \blkdiag\left( \frac{\hat{\ext}_{k|k-1}^{(\ell,1)}}{n_{z,k}^{(\theta,1)}} , \ldots , \frac{\hat{\ext}_{k|k-1}^{(\ell,N_{s,k})}}{n_{z,k}^{(\theta,N_{s,k})}} \right),\\
\hat{\ext}_{k|k-1}^{(\ell,i)} & = \frac{V_{k|k-1}^{(\ell,i)}}{v_{k|k-1}^{(\ell,i)}-2d-2}.
\end{align}
\end{subequations}

\subsubsection{Random matrices}
The corrected parameters are
\begin{subequations}
\begin{align}
v_{k|k}^{(\theta,\ell,i)} = & v_{k|k-1}^{(\ell,i)} + n_{z,k}^{(\theta,i)},\\
V_{k|k}^{(\theta,\ell,i)} = & V_{k|k-1}^{(\ell,i)} + Z_{k}^{(\theta,i)} + N_{k|k-1}^{(\theta,\ell,i)},\\
%\end{align}
%where
%\begin{align}
N_{k|k-1}^{(\theta,\ell,i)} = & \left(\hat{\ext}_{k|k-1}^{(\ell,i)}\right)^{\frac{1}{2}} \left(S_{k}^{(\theta,\ell,i)}\right)^{-\frac{1}{2}} \varepsilon_{k|k-1}^{(\theta,\ell,i)} \nonumber \\
& \times \left(\varepsilon_{k|k-1}^{(\theta,\ell,i)}\right)^{\tp} \left(S_{k}^{(\theta,\ell,i)}\right)^{-\frac{\tp}{2}} \left(\hat{\ext}_{k|k-1}^{(\ell,i)}\right)^{\frac{\tp}{2}},\\
\varepsilon_{k|k-1}^{(\theta,\ell,i)} = & \bar{\sz}_{k}^{(\theta,i)}-H_{k}^{(i)}m_{k|k-1}^{(\ell)},\\
S_{k}^{(\theta,\ell,i)} = & H_{k}^{(i)}P_{k|k-1}^{(\ell)}\left(H_{k}^{(i)}\right)^{\tp} + \frac{\hat{\ext}_{k|k-1}^{(\ell,i)}}{n_{z,k}^{(\theta,i)}},
\end{align}
\end{subequations}
where the matrix square-roots are computed using, \egp, Cholesky factorization.

\subsubsection{Weights}
The weights are computed as
\begin{subequations}
\begin{align}
w_{k|k}^{(\theta,\ell)} = & \frac{w_{k|k-1}^{(\ell)} \prod_{i=1}^{N_{s,k}}\mathcal{L}_{k}^{(\theta,\ell,i)}}{\sum_{\theta'\in\bar{\setTheta}} \sum_{\ell'=1}^{J_{k|k-1}} w_{k|k-1}^{(\ell')} \prod_{i'=1}^{N_{s,k}}\mathcal{L}_{k}^{(\theta',\ell',i')}},\\
%\end{align}
%where the likelihood corresponding to the $i$th subobject under assocation event $\theta$ is
%\begin{align}
\mathcal{L}_{k}^{(\theta,\ell,i)} = & \frac{\Gamma\left(\alpha_{k|k}^{(\theta,\ell,i)}\right)}{\Gamma\left(\alpha_{k|k-1}^{(\ell,i)}\right)} \frac{\left(\beta_{k|k-1}^{(\ell,i)}\right)^{\alpha_{k|k-1}^{(\ell,i)}}}{\left(\beta_{k|k}^{(\theta,\ell,i)}\right)^{\alpha_{k|k}^{(\theta,\ell,i)}}} \nonumber \\
& \times \frac{ \left(n_{z,k}^{(\theta,i)} \pi^{n_{z,k}^{(\theta,i)}} \right)^{-\frac{d}{2}} 2^{-\frac{n_{z,k}^{(\theta,i)} (d-1)}{2}} }{\left|\left(\hat{\ext}_{k|k-1}^{(\ell,i)}\right)^{-\frac{1}{2}}S_{k}^{(\theta,\ell,i)} \left(\hat{\ext}_{k|k-1}^{(\ell,i)}\right)^{-\frac{\tp}{2}} \right|^{\frac{1}{2}}} \nonumber \\
& \times  \frac{\Gamma_{d}\left(\frac{v_{k|k}^{(\theta,\ell,i)}-d-1}{2}\right)}{\Gamma_{d}\left(\frac{v_{k|k-1}^{(\ell,i)}-d-1}{2}\right)} \frac{\left|V_{k|k-1}^{(\ell,i)} \right|^{\frac{v_{k|k-1}^{(\ell,i)}-d-1}{2}}}{ \left|V_{k|k}^{(\theta,\ell,i)} \right|^{\frac{v_{k|k}^{(\theta,\ell,i)}-d-1}{2}} }.
\end{align}
\end{subequations}

\subsubsection{Predicted likelihoods}
The predicted likelihoods are
\begin{subequations}
\begin{align}
& p^{(\ell)}\left(\setZ_{k}|\theta,\setZ^{k-1}\right) = \prod_{i=1}^{N_{s,k}} \mathcal{L}_{k}^{(\theta,\ell,i)} \\
& p\left(\setZ_{k}|\setZ^{k-1}\right) = \frac{1}{|\bar{\setTheta}|} \sum_{\theta\in\bar{\setTheta}} \sum_{\ell=1}^{J_{k|k-1}} w^{(\ell)}_{k|k-1} \prod_{i=1}^{N_{s,k}} \mathcal{L}_{k}^{(\theta,\ell,i)}.
\end{align}
\end{subequations}

\begin{table}[tb]
\caption{Change main subobject}
\label{tab:ChangeMainSubobject}
%\centering
\vspace{-0.5cm}
\rule[0pt]{\columnwidth}{1pt}
\vspace{-0.3cm}
\begin{algorithmic}[1]
\STATE \textbf{Input:} Component $\GGIWpdf{\xi_{k}}{\zeta_{k|k}}$.
\STATE Center: $\hat{\mathbf{p}}_{k|k}^{(c)} = \sum_{i=1}^{N_{s,k}}\hat{\mathbf{p}}_{k|k}^{(i)}$, where $\hat{\mathbf{p}}_{k|k}^{(i)} = \E\left[\left.\mathbf{p}_{k}^{(i)}\right|\setZ^{k}\right]$.
\STATE Distances to center: $\delta_{k|k}^{(i)} = \left\| \hat{\mathbf{p}}_{k|k}^{(c)} - \hat{\mathbf{p}}_{k|k}^{(i)} \right\|_{2}$.
\STATE Closest to center: $j = \argmin{i} \delta_{k|k}^{(i)}$
\IF{ $j \neq 1$ }
	\STATE Change main subobject from $1$ to $j$.
	%\STATE $\left(\tilde{\alpha}_{k|k}^{(1)},\tilde{\beta}_{k|k}^{(1)}\right) = \left({\alpha}_{k|k}^{(j)},{\beta}_{k|k}^{(j)}\right)$, $\left(\tilde{\alpha}_{k|k}^{(j)},\tilde{\beta}_{k|k}^{(j)}\right) = \left({\alpha}_{k|k}^{(1)},{\beta}_{k|k}^{(1)}\right)$, and $\left(\tilde{\alpha}_{k|k}^{(i)},\tilde{\beta}_{k|k}^{(i)}\right) = \left({\alpha}_{k|k}^{(i)},{\beta}_{k|k}^{(i)}\right)$ for $i\neq 1,j$.
	\STATE $\tilde{m}_{k|k} = A_{1,j}m_{k|k}$, $\tilde{P}_{k|k} = A_{1,j}P_{k|k}A_{1,j}^{\tp}$, where $A_{1,j}$ is a permutation matrix that changes place between subobjects $1$ and $j$.
	%\STATE $\left(\tilde{v}_{k|k}^{(1)},\tilde{V}_{k|k}^{(1)}\right) = \left({v}_{k|k}^{(j)},{V}_{k|k}^{(j)}\right)$, $\left(\tilde{v}_{k|k}^{(j)},\tilde{V}_{k|k}^{(j)}\right) = \left({v}_{k|k}^{(1)},{V}_{k|k}^{(1)}\right)$, and $\left(\tilde{v}_{k|k}^{(i)},\tilde{V}_{k|k}^{(i)}\right) = \left({v}_{k|k}^{(i)},{V}_{k|k}^{(i)}\right)$ for $i\neq 1,j$.
	
	\STATE $\left(\tilde{\alpha}_{k|k}^{(1)},\tilde{\beta}_{k|k}^{(1)},\tilde{v}_{k|k}^{(1)},\tilde{V}_{k|k}^{(1)}\right) = \left({\alpha}_{k|k}^{(j)},{\beta}_{k|k}^{(j)},{v}_{k|k}^{(j)},{V}_{k|k}^{(j)}\right)$, 
	
	$\left(\tilde{\alpha}_{k|k}^{(j)},\tilde{\beta}_{k|k}^{(j)},\tilde{v}_{k|k}^{(j)},\tilde{V}_{k|k}^{(j)}\right) = \left({\alpha}_{k|k}^{(1)},{\beta}_{k|k}^{(1)},{v}_{k|k}^{(1)},{V}_{k|k}^{(1)}\right)$, and
	
	$\left(\tilde{\alpha}_{k|k}^{(i)},\tilde{\beta}_{k|k}^{(i)},\tilde{v}_{k|k}^{(i)},\tilde{V}_{k|k}^{(i)}\right) = \left({\alpha}_{k|k}^{(i)},{\beta}_{k|k}^{(i)},{v}_{k|k}^{(i)},{V}_{k|k}^{(i)}\right)$ for $i\neq 1,j$.
	
	
	\STATE $\tilde{\zeta}_{k|k} = \left( \left\{\tilde{\alpha}_{k|k}^{(i)},\tilde{\beta}_{k|k}^{(i)} \right\}_{i=1}^{N_{s,k}},\tilde{m}_{k|k},\tilde{P}_{k|k}, \left\{\tilde{v}_{k|k}^{(i)},\tilde{V}_{k|k}^{(i)} \right\}_{i=1}^{N_{s,k}} \right) $.
\ELSE
	\STATE No change of main subobject, $\tilde{\zeta}_{k|k} = \zeta_{k|k}$
\ENDIF
\STATE \textbf{Output:} Component $\GGIWpdf{\xi_{k}}{\tilde{\zeta}_{k|k}}$.
\end{algorithmic}
\vspace{-0.1cm}
\rule[0pt]{\columnwidth}{1pt}
\end{table}

\subsection{Mixture reduction}
With $J_{k|k}$ components, $M_{k}$ motion models and $|\bar{\setTheta}|$ association events there are $J_{k+1|k+1} = |\bar{\setTheta}| M_{k} J_{k|k}$ components after one iteration of prediction and correction. Mixture reduction is used in each iteration after the correction step to keep the number of components at a tractable level.
%With time the number of components becomes too large and approximations are therefore necessary. In this paper, after the correction step mixture reduction is performed to reduce the number of hypotheses.
Hypotheses with weights lower than a threshold $\tau$ are pruned and the weights are re-normalized. Merging is then performed on the mixture, where we have used a combination of the gamma mixture merging from \cite{GranstromO:2012c} and the Gaussian inverse Wishart merging from \cite{GranstromO:2012e}. Note that when multiple motion models are used, merging is only performed within the same motion modes, and not across the motion modes.

\subsection{Change of main subobject}
When the motion model $f^{(m)}(\sx_k)$ is linear it does not matter which subobject is defined as the main one. However with non-linear motion, \egp a coordinated turn, it is typically most intuitive to define the main subobject as the subobject closest to the center of the overall extended target. Here center is defined as the mean of the subobjects' positions.

It can happen that an estimate hypothesis $p^{(\ell)}(\cdot)$ initialized by the method in Table~\ref{tab:multi_ellipse_initilization}, after a couple of prediction-correction-iterations, converges to a configuration where the main subobject is not the one closest to the center. In this case another subobject can be defined as the main one, an operation that corresponds to a simple re-ordering of the measurement rates and random matrices, and a simple linear transformation of the kinematic vector. Note that any non-linear dynamics could be re-defined too, however empirically we have found that it is sufficient to keep the kinematics $\mathbf{c}_{k}$ constant during the change of main subobject. The algorithm that is used is given in Table~\ref{tab:ChangeMainSubobject}.

An important special case is when there are two subobjects, since in this case both subobjects are equally close to the center. In this case the main subobject is arbitrarily assigned upon initialization, and the change of main subobject is never utilized in the filter recursion.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%
%%%%% Comparison
%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Comparison to alternative models}
In this section we discuss and compare the proposed model to other extended target models that are available in the litterature. Models for specific geometric shapes, \egp sticks, circles, ellipse and rectangles can be found in \cite{GilholmS:2005,BaumFH:2012,Boersetal:06a,PetrovMGA:2011,Koch:2008,BaumNH:2010,GranstromLO:2011,ZhuHL:2011,ReuterD:2011,DegermanWS:2011,LanRL:2012,ReuterWD:2012}. Because these models consider specific shapes we do not compare to them further. Three models capable of handling general and irregular shapes can be found in \cite{LundquistGO:11,BaumH:2011,LanRL:2014}. The model in \cite{LundquistGO:11} considers measurements that are spread along the outline of the target's shape (e.g. laser range measurements), and in this paper we consider measurement that are spread across the target's surface. The model in \cite{LundquistGO:11} is thus not applicable to the scenarios considered here. In the next section we present simulation results that compare the proposed model to the two models presented in \cite{BaumH:2011,LanRL:2014}, and in the remainder of this section we elaborate on the theoretical similarities and differences between the proposed model and the ones from \cite{BaumH:2011,LanRL:2014}.

\subsection{Star-Convex model \cite{BaumH:2011}}
The Star-Convex Random Hypersurface Model \cite{BaumH:2011}, denoted M1, parametrizes the boundary of the target shape as a radial function. The radial function is parameterized using Fourier series, and the Fourier coefficients are estimated. With more coefficients, the radial function has more degrees of freedom and increasingly complex shapes can be described. In a sense, this is analogous to how a larger number of subobjects can describe a more complex shape. Model M1 does not decompose the extended target into subobjects, and thus does not need to solve the measurement-to-subobject association problem.

\subsection{Multi-Ellipse model \cite{LanRL:2014}}
The model in \cite{LanRL:2014}, denoted M2, models the target using multiple elliptic subobjects, see \eqref{eq:LanRongLiStateModel}. The proposed model is very similar to model M2 as they both extend the random matrix framework \cite{Koch:2008,FeldmannFK:2011} to model the subobjects. The two models also have the following differences:

\subsubsection{Measurement rates}
The proposed model includes a model of the number of detections per subobject per timestep and estimates the measurement rates for each subobject, which M2 does not.

\subsubsection{Main subobject}
The proposed model defines one of the subobjects as the main subobject, around which remaining subobjects are located. %With linear motion models it does not matter which subobject is the main one. However, with non-linear motion models this becomes important, \egp when the target performs a rotating maneuver, since the position of the main subobject is the point around which the target rotates. 
In comparison, M2 does not define one of the subobjects as the main one.

\subsubsection{Position covariance and unified kinematics}
M2 is based on the random matrix model \eqref{eq:Koch_random_matrix_model}, see \eqref{eq:LanRongLiStateModel}, while the proposed model is based on the random matrix model \eqref{eq:Feldmann_random_matrix_model}. This difference is fundamental, because it is what allows the positions and kinematics of all subobjects to the modeled as a single random vector. Due to the form of the Gaussian covariances in \eqref{eq:Koch_random_matrix_model} ($P\otimes X$), under this model the state vectors of the subobjects cannot easily be treated as a single random vector. Modeling with a single random vector improves the overall modeling in the following ways:
\begin{enumerate}
	\item The proposed model estimates unified kinematics $\mathbf{c}_{k}$ (\iep a single velocity, a single turn-rate etc for the extended target as a whole) for all subobjects', and if necessary individual subobject kinematics can be included in $\mathbf{c}_{k}$. In comparison, M2 estimates individual velocities and accelerations for the subobjects. A multiple motion model framework is used in M2 where there are some common kinematics via the process noise parameters, however unified kinematics are not estimated.
	\item The proposed model maintains a full covariance matrix for the subobjects' positions and the kinematics, i.e. the dependencies between the subobjects' positions and the kinematics are modeled and estimated. In comparison, M2 does not model the dependencies between the subobjects' kinematic states $\sx_{k}^{(i)}$, cf \eqref{eq:LanRongLiStateModel}.
\end{enumerate}

\begin{figure*}[htbp]
\centerline{
\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_true_tracks_xy} \label{fig:multi_ellipse_true_tracks_xy}}
\hfill
\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_true_tracks_v} \label{fig:multi_ellipse_true_tracks_v}}
\hfill
\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_true_tracks_omega} \label{fig:multi_ellipse_true_tracks_omega}}
}
\caption{True target trajectory, initial position is origin. Left: $x,y$-position. Middle: velocity. Right: turn-rate.}
\label{fig:multi_ellipse_true_tracks}
\end{figure*}

\subsubsection{Kinematics-extension-dependence}
The proposed model does not model the dependence between the position and the extension of a subobject, which M2 does. Further, under an assumed linear Gaussian measurement model the proposed model's measurement update is approximate, while M2's is exact. Both these differences follow from the proposed model being based on the random matrix model \eqref{eq:Feldmann_random_matrix_model} and M2 being based on the random matrix model \eqref{eq:Koch_random_matrix_model}.

With regards to the dependence not being modeled, the measurement update and the prediction update provides for the necessary interdependence between the kinematic state estimate and the extension state estimate. This is analogous to how the measurement update \cite{FeldmannFK:2011} and the prediction update \cite{GranstromO:2014} provide for interdependence in case the target is modeled using a single ellipse. In simulation comparisons for single elliptic targets, models based on \cite{FeldmannFK:2011} outperform models based on \cite{Koch:2008}, see \cite{GranstromO:2014}.

We conclude the comparison by noting that different models should not be judged and compared only on their theoretical properties, but also on their practical properties. In the next section we present a simulation study that compares the practical performance of the proposed model and models M1 and M2.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%
%%%%% Simulation results
%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Simulation results}
\label{sec:simulation_results}
%In this section we present results from a simulation study. The presented extended target tracking method is compared to a version of the method presented in \cite{LanRL:2014}, where measurement rate estimation has been added as outlined in \cite{GranstromO:2012c}.

\subsection{Target extraction and performance evaluation}
\label{sec:performance_eval}
To extract a target estimate from a mixture \eqref{eq:mixtureGGIW_extended_target_state}, merging is first performed, this time across the motion modes. Expected values of the measurement rates, positions and extension matrices are then computed w.r.t. the component with the highest weight $w_{k|k}^{(\ell)}$. Both the predicted estimate $\hat{\xi}_{k|k-1}$ and the filtered estimate $\hat{\xi}_{k|k}$ is compared to the true target state $\xi_{k}$. The following error metrics are used for the measurement rates, subobject positions, and random matrices,
\begin{subequations}
\begin{align}
d^{\gamma}_{k|k} = & \sum_{i=1}^{N_{s,k}} \left| \gamma_{k}^{(i)} - \hat{\gamma}_{k|k}^{(\pi(i))} \right|, \ \hat{\gamma}_{k|k}^{(i)}=\E\left[\left.\gamma_{k}^{(i)}\right|\setZ^{k}\right]\\
d^{\mathbf{p}}_{k|k} = & \sum_{i=1}^{N_{s,k}} \left\| \mathbf{p}_{k}^{(i)} - \hat{\mathbf{p}}_{k|k}^{(\pi(i))}\right\|_{2}, \ \hat{\mathbf{p}}_{k|k}^{(i)} = \E\left[ \left. \mathbf{p}_{k}^{(i)} \right| \setZ^{k} \right] \\
d^{\ext}_{k|k} = & \sum_{i=1}^{N_{s,k}} \left\| \ext_{k}^{(i)} - \hat{\ext}_{k|k}^{(\pi(i))} \right\|_{F}, \hat{\ext}_{k|k}^{(i)}=\E\left[\left.\ext_{k}^{(i)}\right|\setZ^{k}\right]
\end{align}
\end{subequations}
A subobject-to-subobject association $\pi(i)$ is obtained by minimizing $d^{\mathbf{p}}_{k|k}$. Because $\gamma_{k}^{(i)}$, $\mathbf{p}_{k}^{(i)}$ and $\ext_{k}^{(i)}$ all have different units we refrain from computing an overall metric for the extended target state $\xi_{k}$.

\subsection{True tracks and setup}
We simulate both stationary and moving targets. In a comparison of shape estimation, stationary targets are used because we wish to emphasize the shape estimation, not the motion estimation. For the moving target, the target trajectory that was simulated is shown in \figurename~\ref{fig:multi_ellipse_true_tracks}; the true position is shown in \figurename~\ref{fig:multi_ellipse_true_tracks_xy} and the corresponding speed and turn-rate is shown in \figurename~\ref{fig:multi_ellipse_true_tracks_v} and \figurename~\ref{fig:multi_ellipse_true_tracks_omega}.

Three different $d=2$ dimensional extended target shapes were simulated. The first is T-shaped, with measurements generated by uniformly sampling across the shape and adding Gaussian noise with covariance $\mathbf{R}=\mathbf{I}_{2}$. The measurement rate was $4\gamma_{0}$. The other two shapes consist of three and two subobjects, respectively. The shape of the targets are consistent with the examples given in \figurename~\ref{fig:multiple_ellipse_example_target}, \iep the shape resembles that of an airplane and of the letter V, respectively. For these two shapes, the measurement model \eqref{eq:Linear_Gaussian_Measurements} was used. For the plane-like target, for the subobject that corresponds to the fuselage the measurement rate was $2\gamma_{0}$, and the extension matrix was $\ext=\diag{[10^2 \ , \ 2^2]}$. For the subobjects that correspond to the wings the measurement rates were $\gamma_{0}$, and the extension matrices were $\ext=\diag{[5^2 \ , \ 1^2]}$. For the V-shaped target, the subobjects both had measurement rates $\gamma_{0}$ and extension matrices $\ext=\diag{[20^2\ , \ 1^2]}$. The scenarios were simulated for different values of $\gamma_{0}$: $2$, $5$ and $20$.

For the presented filter, a constant turn-rate (\ctmod) motion model $f_{\mathbf{p},\mathbf{c}}^{(m)}(\cdot)$ with polar velocity, see \cite[Eq. 75]{RongLiJ:2003}, was used for the overall position and kinematics. In this case the unified kinematics are given by $\mathbf{c}_{k} = \begin{bmatrix} \cv_{k} & \phi_{k} & \omega_{k} \end{bmatrix}^{\tp}$, where $\cv_{k}$ is the speed, $\phi_{k}$ is the heading and $\omega_{k}$ is the turn-rate. For this type of motion model the time evolution of the subobject offsets is
\begin{align}
\mathbf{d}_{k+1}^{(i)} = f_{\mathbf{d}}^{(m)}\left(\mathbf{d}_{k}^{(i)},\mathbf{c}_{k}\right) = R \left(T\omega_{k}\right) \mathbf{d}_{k}^{(i)}
\end{align}
where $T$ is the sampling time and $R \left(\cdot\right)$ is a rotation matrix. The matrix transformation function is also a rotation matrix, $M^{(m)}\left(\sx_k\right) = R \left(T\omega_{k}\right)$.
%\begin{align}
%M\left(\sx_k\right) = R \left(T\omega_{k}\right).
%\end{align}
Two \ctmod models were implemented, one with small process noise corresponding to non-maneuver, and one with larger process noise corresponding to maneuver. The transition probabilities were set to $95\%$ probability to stay in the same mode, and $5\%$ probability for mode switch.
%The transition probability matrix was set to
%\begin{align}
%\begin{bmatrix} 0.95 & 0.05 \\ 0.05 & 0.95 \end{bmatrix}.
%\end{align} 

\begin{figure*}[!t]
\centerline{
\subfloat{\includegraphics[width=0.125\textwidth]{multi_ellipse_shape_comparison_gam_2_ex_1} \label{fig:multi_ellipse_shape_comparison_gam_2_ex_1}}
\hfill
\subfloat{\includegraphics[width=0.125\textwidth]{multi_ellipse_shape_comparison_gam_2_ex_2} \label{fig:multi_ellipse_shape_comparison_gam_2_ex_2}}
\hfill
\subfloat{\includegraphics[width=0.125\textwidth]{multi_ellipse_shape_comparison_gam_2_ex_3} \label{fig:multi_ellipse_shape_comparison_gam_2_ex_3}}
\hfill
\subfloat{\includegraphics[width=0.125\textwidth]{multi_ellipse_shape_comparison_gam_2_ex_4} \label{fig:multi_ellipse_shape_comparison_gam_2_ex_4}}
\hfill
\subfloat{\includegraphics[width=0.125\textwidth]{multi_ellipse_shape_comparison_Tshape_gam_2_ex_1} \label{fig:multi_ellipse_shape_comparison_Tshape_gam_2_ex_1}}
\hfill
\subfloat{\includegraphics[width=0.125\textwidth]{multi_ellipse_shape_comparison_Tshape_gam_2_ex_2} \label{fig:multi_ellipse_shape_comparison_Tshape_gam_2_ex_2}}
\hfill
\subfloat{\includegraphics[width=0.125\textwidth]{multi_ellipse_shape_comparison_Tshape_gam_2_ex_3} \label{fig:multi_ellipse_shape_comparison_Tshape_gam_2_ex_3}}
\hfill
\subfloat{\includegraphics[width=0.125\textwidth]{multi_ellipse_shape_comparison_Tshape_gam_2_ex_4} \label{fig:multi_ellipse_shape_comparison_Tshape_gam_2_ex_4}}
}
\centerline{
\subfloat{\includegraphics[width=0.125\textwidth]{multi_ellipse_shape_comparison_gam_20_ex_1} \label{fig:multi_ellipse_shape_comparison_gam_20_ex_1}}
\hfill
\subfloat{\includegraphics[width=0.125\textwidth]{multi_ellipse_shape_comparison_gam_20_ex_2} \label{fig:multi_ellipse_shape_comparison_gam_20_ex_2}}
\hfill
\subfloat{\includegraphics[width=0.125\textwidth]{multi_ellipse_shape_comparison_gam_20_ex_3} \label{fig:multi_ellipse_shape_comparison_gam_20_ex_3}}
\hfill
\subfloat{\includegraphics[width=0.125\textwidth]{multi_ellipse_shape_comparison_gam_20_ex_4} \label{fig:multi_ellipse_shape_comparison_gam_20_ex_4}}
\hfill
\subfloat{\includegraphics[width=0.125\textwidth]{multi_ellipse_shape_comparison_Tshape_gam_20_ex_1} \label{fig:multi_ellipse_shape_comparison_Tshape_gam_20_ex_1}}
\hfill
\subfloat{\includegraphics[width=0.125\textwidth]{multi_ellipse_shape_comparison_Tshape_gam_20_ex_2} \label{fig:multi_ellipse_shape_comparison_Tshape_gam_20_ex_2}}
\hfill
\subfloat{\includegraphics[width=0.125\textwidth]{multi_ellipse_shape_comparison_Tshape_gam_20_ex_3} \label{fig:multi_ellipse_shape_comparison_Tshape_gam_20_ex_3}}
\hfill
\subfloat{\includegraphics[width=0.125\textwidth]{multi_ellipse_shape_comparison_Tshape_gam_20_ex_4} \label{fig:multi_ellipse_shape_comparison_Tshape_gam_20_ex_4}}
}
\caption{Example results for tracking of a plane-shaped and a T-shaped stationary target. Top row $\gamma_{0}=2$, bottom row $\gamma_{0}=20$. True target (gray area) compared to proposed model (solid blue line, main subobject indicated by thicker line), model M2 (dashed orange line), model M1 (dash-dotted red line). Time step $k$ written in lower left corner. As expected convergence is much faster when there are more measurements (i.e. higher $\gamma_{0}$).}
\label{fig:multi_ellipse_shape_comparison}
\end{figure*}

\begin{figure*}[!t]
\centerline{
\subfloat[$\gamma_{0}=2$]{\includegraphics[width=0.33\textwidth]{multi_ellipse_planeShape_MC_boxplot_gam_2} \label{fig:multi_ellipse_planeShape_MC_boxplot_gam_2}}
\vline
\subfloat[$\gamma_{0}=5$]{\includegraphics[width=0.33\textwidth]{multi_ellipse_planeShape_MC_boxplot_gam_5} \label{fig:multi_ellipse_planeShape_MC_boxplot_gam_5}}
\vline
\subfloat[$\gamma_{0}=20$]{\includegraphics[width=0.33\textwidth]{multi_ellipse_planeShape_MC_boxplot_gam_20} \label{fig:multi_ellipse_planeShape_MC_boxplot_gam_20}}
}
\caption{Estimation errors for plane-shaped target. Proposed model in blue, model M2 in orange. x-labels F and P denote filter errors $d_{k|k}$ and prediction errors $d_{k|k-1}$. On each box, central mark is median, edges of box are 25th and 75th percentiles, whiskers extend to most extreme datapoints the algorithm considers to be not outliers.}
\label{fig:multi_ellipse_planeShape_MC_boxplot}
\end{figure*}

The filter parameters that were used in the implementation are listed in Table~\ref{tab:FilterParameters}.

\begin{table}[htbp]
	\centering
	\caption{Parameters for proposed method}
		\begin{tabular}{lcc}
			Parameter &  & Value \\
			\hline
			Sample time & $T$ & 1 \\
			Number of initial hypotheses & $N_p$ & $2(N_{s,k}-1)$ \\
			Initial kinematics & $\mathbf{c}_{0}$ & $\mathbf{0}_{3\times 1}$ \\
			Initial covariance & $P_{0}$ & $10^{2}\mathbf{I}_{n_x}$ \\
			Measurement rate initial mean & $e$ & 15 \\
			Measurement rate initial variance & $v$ & 10 \\
			%Number of motion models & $M_{k}$ & $2$ \\
			Measurement rate prediction factor & $\eta_{k}^{(m)}$ & $1.05$, $\forall m$ \\
			%Motion model & $f^{(m)}(\sx_{k})$ & ? \\
			%Matrix transformation & $M^{(m)}(\sx_{k})$ & ? \\
			Prediction degrees of freedom & $n_{k+1}^{(m)}$ & $100$, $\forall m$ \\
			Pruning threshold & $\tau$ & $0.01$
		\end{tabular}
	\label{tab:FilterParameters}
\end{table}

%The proposed model is compared to the Star-Convex Random Hypersurface Model (\textsc{scrhm}) \cite{BaumH:2011} and the multi ellipse model from \cite{LanRL:2014}. For the \textsc{scrhm} an implementation available online was used\footnote{Thanks to M. Baum and R. Sandkuehler for providing code.\\ {\tiny\texttt{http://www.cloudrunner.eu/algorithm/12/random-hypersurface-model/version/2/}}}.

For M1 an implementation available online was used\footnote{Thanks to M. Baum and R. Sandkuehler for providing code.\\ {\tiny\texttt{http://www.cloudrunner.eu/algorithm/12/random-hypersurface-model/version/2/}}}. Model M2 was implemented as instructed in Section VI ``Simulation studies'' in \cite{LanRL:2014}, and augmented to include estimation of the measurement rates, see \cite{GranstromO:2012c} for details. The method was parametrized with three motion models: the first corresponds to constant velocity motion; the second corresponds to a turn with turn-rate $\psi$; the third corresponds to a turn with turn-rate $-\psi$. In the simulations the parameter $\psi$ was set to $5$ degrees per second. For the lower measurement rates ($\gamma_{0}=2$ and $\gamma_{0}=5$) sometimes during the maneuvers one of the subobjects would separate from the other two subobjects. To correct this, any hypothesis $p^{(\ell)}(\cdot)$ with a subobject more than $50$ meters from the other subobjects was deleted. In case all hypotheses were deleted, the target estimate was reinitialized.

%\begin{figure}[!t]
%\centerline{
%\subfloat{\includegraphics[width=0.25\columnwidth]{multi_ellipse_shape_comparison_Tshape_gam_2_ex_1} \label{fig:multi_ellipse_shape_comparison_Tshape_gam_2_ex_1}}
%\hfill
%\subfloat{\includegraphics[width=0.25\columnwidth]{multi_ellipse_shape_comparison_Tshape_gam_2_ex_2} \label{fig:multi_ellipse_shape_comparison_Tshape_gam_2_ex_2}}
%\hfill
%\subfloat{\includegraphics[width=0.25\columnwidth]{multi_ellipse_shape_comparison_Tshape_gam_2_ex_3} \label{fig:multi_ellipse_shape_comparison_Tshape_gam_2_ex_3}}
%\hfill
%\subfloat{\includegraphics[width=0.25\columnwidth]{multi_ellipse_shape_comparison_Tshape_gam_2_ex_4} \label{fig:multi_ellipse_shape_comparison_Tshape_gam_2_ex_4}}
%}
%\centerline{
%\subfloat{\includegraphics[width=0.25\columnwidth]{multi_ellipse_shape_comparison_Tshape_gam_20_ex_1} \label{fig:multi_ellipse_shape_comparison_Tshape_gam_20_ex_1}}
%\hfill
%\subfloat{\includegraphics[width=0.25\columnwidth]{multi_ellipse_shape_comparison_Tshape_gam_20_ex_2} \label{fig:multi_ellipse_shape_comparison_Tshape_gam_20_ex_2}}
%\hfill
%\subfloat{\includegraphics[width=0.25\columnwidth]{multi_ellipse_shape_comparison_Tshape_gam_20_ex_3} \label{fig:multi_ellipse_shape_comparison_Tshape_gam_20_ex_3}}
%\hfill
%\subfloat{\includegraphics[width=0.25\columnwidth]{multi_ellipse_shape_comparison_Tshape_gam_20_ex_4} \label{fig:multi_ellipse_shape_comparison_Tshape_gam_20_ex_4}}
%}
%\caption{Example results for tracking of a T-shaped stationary target. Top row $\gamma_{0}=2$, bottom row $\gamma_{0}=20$. True target (gray area) compared to proposed model (solid blue line, main subobject indicated by thicker line), model from \cite{LanRL:2014} (dashed orange line), model from \cite{BaumH:2011} (dash-dotted red line). Time step $k$ written in lower left corner. As expected convergence is much faster when there are more measurements (i.e. higher $\gamma_{0}$).}
%\label{fig:multi_ellipse_shape_comparison_Tshape}
%\end{figure}

\begin{figure*}[htp]
\centerline{
\subfloat[$\gamma_{0}=2$]{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_boxplot_gam_2} \label{fig:multi_ellipse_Vshape_MC_boxplot_gam_2}}
\vline
\subfloat[$\gamma_{0}=5$]{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_boxplot_gam_5} \label{fig:multi_ellipse_Vshape_MC_boxplot_gam_5}}
\vline
\subfloat[$\gamma_{0}=20$]{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_boxplot_gam_20} \label{fig:multi_ellipse_Vshape_MC_boxplot_gam_20}}
}
\caption{Estimation errors for V-shaped target. Proposed model in blue, model M2 in orange. x-labels F and P denote filter errors $d_{k|k}$ and prediction errors $d_{k|k-1}$. On each box, central mark is median, edges of box are 25th and 75th percentiles, whiskers extend to most extreme datapoints the algorithm considers to be not outliers.}
\label{fig:multi_ellipse_Vshape_MC_boxplot}
\end{figure*}

\subsection{Results: stationary target}
The T-shaped and the plane-shaped targets were simulated for measurement rates $\gamma_{0}=2$ and $\gamma_{0}=20$. A comparison of the proposed model and the models M1 and M2 is shown in \figurename~\ref{fig:multi_ellipse_shape_comparison}. As expected all three methods converge much faster when there are more measurements (\iep higher $\gamma_{0}$). All three methods give reasonable results, however the proposed model and the M2 are closer to the ground truth than M1. Because of this, for the moving target we only compare the proposed model and model M2.

\subsection{Results: moving target}
The plane-shaped and the V-shaped targets were simulated for $\gamma_{0}=2$, $\gamma_{0}=5$, and $\gamma_{0}=20$. For each value of the measurement rate $\gamma_{0}$ the scenarios were simulated $10^{3}$ times.
For the plane-shaped target (three subobjects) the filter errors $d_{k|k}$ and prediction errors $d_{k|k-1}$ are shown in \figurename~\ref{fig:multi_ellipse_planeShape_MC_boxplot}, for the V-shaped target (two subobjects) the results are shown in \figurename~\ref{fig:multi_ellipse_Vshape_MC_boxplot}.
%For the plane-shaped target (three subobjects) the filter errors $d_{k|k}$ and prediction errors $d_{k|k-1}$ are shown in \figurename~\ref{fig:multi_ellipse_MC_gam_filter} and \figurename~\ref{fig:multi_ellipse_MC_gam_pred}, for the V-shaped target (two subobjects) the results are shown in \figurename~\ref{fig:multi_ellipse_Vshape_MC_gam_filter} and \figurename~\ref{fig:multi_ellipse_Vshape_MC_gam_pred}.
Example filter and prediction outputs for the plane-shaped target for $\gamma_{0}=5$ are shown in \figurename~\ref{fig:multi_ellipse_Example_Results}. From the results the following observations can be made:
\begin{itemize}
	\item Both the prediction errors and the filter errors are smaller for the presented method for all $\gamma_{0}$.
	\item The biggest difference is for the subobject position errors, especially during maneuvers. Note that, even if the estimated random matrices have the correct size and orientation, the subobject positions are more important for the overall extended target extension estimate. The larger the subobject position errors are, the more distorted the overall shape becomes, which can be seen in \figurename~\ref{fig:multi_ellipse_Example_Results}.
\end{itemize}
The lower errors for the presented method, especially the lower position errors, are a direct effect of 
\begin{inparaenum}[\itshape a\upshape)]
	\item using a single state vector for the subobject positions and the kinematics including a full covariance matrix; and 
	\item having unified kinematics for the subobject positions.
\end{inparaenum}

%\begin{figure*}[!t]
%\centerline{
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_MC_gam_2_RMSE_g} \label{fig:multi_ellipse_MC_gam_2_RMSE_g}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_MC_gam_2_RMSE_x} \label{fig:multi_ellipse_MC_gam_2_RMSE_x}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_MC_gam_2_RMSE_ext} \label{fig:multi_ellipse_MC_gam_2_RMSE_ext}}
%}
%\centerline{
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_MC_gam_5_RMSE_g} \label{fig:multi_ellipse_MC_gam_5_RMSE_g}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_MC_gam_5_RMSE_x} \label{fig:multi_ellipse_MC_gam_5_RMSE_x}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_MC_gam_5_RMSE_ext} \label{fig:multi_ellipse_MC_gam_5_RMSE_ext}}
%}
%\centerline{
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_MC_gam_20_RMSE_g} \label{fig:multi_ellipse_MC_gam_20_RMSE_g}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_MC_gam_20_RMSE_x} \label{fig:multi_ellipse_MC_gam_20_RMSE_x}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_MC_gam_20_RMSE_ext} \label{fig:multi_ellipse_MC_gam_20_RMSE_ext}}
%}
%\caption{Filter results for plane-shaped target, mean error (solid lines) $\pm$ one standard deviation (colored areas). Top row $\gamma_{0} = 2$, middle row $\gamma_{0} = 5$, bottom row $\gamma_{0} = 20$. Left, middle and right columns show errors for measurement rates, kinematic state, and random matrices, respectively. Presented method in blue and method from \cite{LanRL:2014} in orange.}
%\label{fig:multi_ellipse_MC_gam_filter}
%\end{figure*}

%\begin{figure*}[!t]
%\centerline{
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_MC_gam_2_RMSE_g_pred} \label{fig:multi_ellipse_MC_gam_2_RMSE_g_pred}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_MC_gam_2_RMSE_x_pred} \label{fig:multi_ellipse_MC_gam_2_RMSE_x_pred}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_MC_gam_2_RMSE_ext_pred} \label{fig:multi_ellipse_MC_gam_2_RMSE_ext_pred}}
%}
%\centerline{
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_MC_gam_5_RMSE_g_pred} \label{fig:multi_ellipse_MC_gam_5_RMSE_g_pred}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_MC_gam_5_RMSE_x_pred} \label{fig:multi_ellipse_MC_gam_5_RMSE_x_pred}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_MC_gam_5_RMSE_ext_pred} \label{fig:multi_ellipse_MC_gam_5_RMSE_ext_pred}}
%}
%\centerline{
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_MC_gam_20_RMSE_g_pred} \label{fig:multi_ellipse_MC_gam_20_RMSE_g_pred}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_MC_gam_20_RMSE_x_pred} \label{fig:multi_ellipse_MC_gam_20_RMSE_x_pred}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_MC_gam_20_RMSE_ext_pred} \label{fig:multi_ellipse_MC_gam_20_RMSE_ext_pred}}
%}
%\caption{Prediction results for plane-shaped target, mean error (solid lines) $\pm$ one standard deviation (colored areas). Top row $\gamma_{0} = 2$, middle row $\gamma_{0} = 5$, bottom row $\gamma_{0} = 20$. Left, middle and right columns show errors for measurement rates, kinematic state, and random matrices, respectively. Presented method in blue and method from \cite{LanRL:2014} in orange.}
%\label{fig:multi_ellipse_MC_gam_pred}
%\end{figure*}

\subsection{Computational complexity}
\label{sec:Computational_Complexity}
The code used in this work was implemented in \matlab and run on a $2.83$GHz Intel Core2 Quad CPU with $3.48$GB of RAM running Windows. Note that the code has not been optimized for speed.

%\footnote{Implementation will be made available online if paper is accepted.}

In each time step approximately $15$ to $25$ different partitions of the set of measurements were computed. The average number of measurement-to-subobject association events are given in Table~\ref{tab:AssociationEvents}. A comparison to the number of association events if there are $\E[n_{z,k}]$ measurements shows that the set of assocation events is reduced by several orders of magnitude. It is noteworthy that for $\gamma_{0}=20$ the reduction in number of association events is by far greatest, yet the estimation errors are smaller for $\gamma_{0}=20$ than for $\gamma_{0}=2$ and $\gamma_{0}=5$. 

\begin{table}[htbp]
	\centering
	\caption{Number of association events, mean $\pm$ standard deviation}
		\begin{tabular}{ccc}
			$\gamma_{0}$ & $\left|\bar{\setTheta}\right|$ & $(N_{s,k})^{\E[n_{z,k}]}$ \\
			\hline
			$2$ & $114 \pm 29$ & $6.6 \times 10^3$ \\
			$5$ & $128 \pm 24$ & $3.5\times 10^{9}$ \\
			$20$ & $130 \pm 23$ & $1.5\times 10^{38}$
		\end{tabular}
	\label{tab:AssociationEvents}
\end{table}

The number of mixture components increase in each time step. However, in the mixture reduction step many components can be pruned, and the remaining components can be merged such that typically only $2$ to $6$ components remain.

%The number of mixture components increase after prediction and correction, and must be reduced in each time step. The average number mixture components are given in Table~\ref{tab:NumberOfMixtureComponents}. The number of components decrease with increasing $\gamma_{0}$ because with more measurements the scenario is less ambiguous. Noteworthy is that the number of components can be reduced drastically in the reduction step without sacrificing performance.
%
%\begin{table}[htpb]
	%\centering
	%\caption{Number components in mixture, mean $\pm$ standard deviation}
		%\begin{tabular}{cccc}
			%$\gamma_{0}$ & After prediction & After correction & After reduction \\
			%\hline
			%$2$ & $6 \pm 3$ & $710 \pm 388$ & $3 \pm 1$ \\
			%$5$ & $4 \pm 2$ & $560 \pm 280$ & $2 \pm 1$ \\
			%$20$ & $3 \pm 2$ & $388 \pm 216$ & $1 \pm 1$
		%\end{tabular}
	%\label{tab:NumberOfMixtureComponents}
%\end{table}

The average cycle times are given in Table~\ref{tab:CycleTimes}. We see that the times for prediction and reduction are independent of $\gamma_{0}$. The time to compute clusters for data association increases when $\gamma_{0}$ increases, because with more measurements it takes more time to cluster them. The correction time decreases when $\gamma_{0}$ increases, because with more measurements the scenario is less ambiguous and the probability mixture typically has fewer components. Note that these two increases/decreases in time offset each other such that the average total cycle time is about $1.7$ seconds for all values of $\gamma_{0}$ that were tested.

The average cycle time for the method from \cite{LanRL:2014} is $2.1 \pm 1.4$ seconds for $\gamma_{0}=2$, $1.7\pm 0.7$ seconds for $\gamma_{0}=5$, and $1.6\pm 0.5$ seconds for $\gamma_{0}=20$. When comparing the cycle times, remember that neither implementation was optimized for speed.

\begin{table}[htbp]
	\centering
	\caption{Cycle times [seconds], mean $\pm$ standard deviation}
		\begin{tabular}{cccccc}
			$\gamma_{0}$ & Prediction & Clusters & Correction & Reduction & Total \\
			\hline
			$2$ & $0.4 \pm 0.2$ & $0.7 \pm 0.1$ & $0.6 \pm 0.4$ & $0.02 \pm 0.02$ & $1.7 \pm 0.6$ \\
			$5$ & $0.3 \pm 0.1$ & $0.9 \pm 0.1$ & $0.5 \pm 0.3$ & $0.01 \pm 0.02$ & $1.8 \pm 0.4$ \\
			$20$ & $0.2 \pm 0.1$ & $1.2 \pm 0.1$ & $0.3 \pm 0.2$ & $0.01 \pm 0.02$ & $1.7 \pm 0.3$
		\end{tabular}
	\label{tab:CycleTimes}
\end{table}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%
%%%%% Conclusions and future work
%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Conclusions and future work}
\label{sec:conclusion_future_work}

The paper has presented an extended target model in which the target extension is modeled using a collection of elliptical subobjects. %A derivation of prediction, correction and predicted likelihood was given, and an implementation was suggested. The proposed model was compared to previous extended target models, specifically the models from \cite{BaumH:2011,LanRL:2014}, both in terms of theoretical properties but also in terms of practical extended target estimation performance using simulated data. 
The simulation results show that the proposed model outperforms previous work.


%\begin{figure*}[!t]
%\centerline{
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_gam_2_RMSE_g} \label{fig:multi_ellipse_Vshape_MC_gam_2_RMSE_g}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_gam_2_RMSE_x} \label{fig:multi_ellipse_Vshape_MC_gam_2_RMSE_x}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_gam_2_RMSE_ext} \label{fig:multi_ellipse_Vshape_MC_gam_2_RMSE_ext}}
%}
%\centerline{
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_gam_5_RMSE_g} \label{fig:multi_ellipse_Vshape_MC_gam_5_RMSE_g}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_gam_5_RMSE_x} \label{fig:multi_ellipse_Vshape_MC_gam_5_RMSE_x}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_gam_5_RMSE_ext} \label{fig:multi_ellipse_Vshape_MC_gam_5_RMSE_ext}}
%}
%\centerline{
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_gam_20_RMSE_g} \label{fig:multi_ellipse_Vshape_MC_gam_20_RMSE_g}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_gam_20_RMSE_x} \label{fig:multi_ellipse_Vshape_MC_gam_20_RMSE_x}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_gam_20_RMSE_ext} \label{fig:multi_ellipse_Vshape_MC_gam_20_RMSE_ext}}
%}
%\caption{Filter results for V-shaped target, mean error (solid lines) $\pm$ one standard deviation (colored areas). Top row $\gamma_{0} = 2$, middle row $\gamma_{0} = 5$, bottom row $\gamma_{0} = 20$. Left, middle and right columns show errors for measurement rates, kinematic state, and random matrices, respectively. Presented method in blue and method from \cite{LanRL:2014} in orange.}
%\label{fig:multi_ellipse_Vshape_MC_gam_filter}
%\end{figure*}

The presented model can be reduced to the cases where either the measurement rates $\gamma_{k}^{(i)}$, the extension matrices $\ext_{k}^{(i)}$, or both, are known. The simulation study considered extended targets, however the model is applicable also to group targets. In case the targets in the group are moving relative to each other, in addition to the unified group movement, individual kinematics can be estimated along with the unified kinematics.

An important topic for future work is to include estimation of the number of subobjects. Interesting lines for future work also include integrating the proposed model into a multiple target algorithm, such as an extended target \phd- or \cphd filter, see \egp \cite{mahler_FUSION_2009_extTarg,LundquistGO:2013}. Furthermore alternative measurement models can be considered, see \egp \cite{FeldmannFK:2011,Orguner:2012,LanRL:2012} for recent work on this topic within the random matrix framework.

%\begin{figure*}[!t]
%\centerline{
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_gam_2_RMSE_g_pred} \label{fig:multi_ellipse_Vshape_MC_gam_2_RMSE_g_pred}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_gam_2_RMSE_x_pred} \label{fig:multi_ellipse_Vshape_MC_gam_2_RMSE_x_pred}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_gam_2_RMSE_ext_pred} \label{fig:multi_ellipse_Vshape_MC_gam_2_RMSE_ext_pred}}
%}
%\centerline{
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_gam_5_RMSE_g_pred} \label{fig:multi_ellipse_Vshape_MC_gam_5_RMSE_g_pred}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_gam_5_RMSE_x_pred} \label{fig:multi_ellipse_Vshape_MC_gam_5_RMSE_x_pred}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_gam_5_RMSE_ext_pred} \label{fig:multi_ellipse_Vshape_MC_gam_5_RMSE_ext_pred}}
%}
%\centerline{
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_gam_20_RMSE_g_pred} \label{fig:multi_ellipse_Vshape_MC_gam_20_RMSE_g_pred}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_gam_20_RMSE_x_pred} \label{fig:multi_ellipse_Vshape_MC_gam_20_RMSE_x_pred}}
%\hfil
%\subfloat{\includegraphics[width=0.33\textwidth]{multi_ellipse_Vshape_MC_gam_20_RMSE_ext_pred} \label{fig:multi_ellipse_Vshape_MC_gam_20_RMSE_ext_pred}}
%}
%\caption{Prediction results for V-shaped target, mean error (solid lines) $\pm$ one standard deviation (colored areas). Top row $\gamma_{0} = 2$, middle row $\gamma_{0} = 5$, bottom row $\gamma_{0} = 20$. Left, middle and right columns show errors for measurement rates, kinematic state, and random matrices, respectively. Presented method in blue and method from \cite{LanRL:2014} in orange.}
%\label{fig:multi_ellipse_Vshape_MC_gam_pred}
%\end{figure*}

%\appendices\section{Derivation of measurement update}
%\appendix
%\label{app:measurement_update}


\begin{figure}[!t]
\centerline{
\subfloat{\includegraphics[width=0.90\columnwidth]{multi_ellipse_Example_Results_Straight_Filter} \label{fig:multi_ellipse_Example_Results_Straight_Filter}}
}
\centerline{
\subfloat{\includegraphics[width=0.90\columnwidth]{multi_ellipse_Example_Results_Filter} \label{fig:multi_ellipse_Example_Results_Filter}}
}
\centerline{
\subfloat{\includegraphics[width=0.90\columnwidth]{multi_ellipse_Example_Results_Straight_Prediction} \label{fig:multi_ellipse_Example_Results_Straight_Prediction}}
}
\centerline{
\subfloat{\includegraphics[width=0.90\columnwidth]{multi_ellipse_Example_Results_Prediction} \label{fig:multi_ellipse_Example_Results_Prediction}}
}
\caption{Example results for $\gamma_{0}=5$ for the trajectory in \figurename~\ref{fig:multi_ellipse_true_tracks}. Top: filtered estimates. Bottom: predicted estimates. Ground truth (gray area), compared to model M2 (dashed orange line) and proposed method (solid blue line, main subobject drawn with thicker line).}
\label{fig:multi_ellipse_Example_Results}
\end{figure}

%\section*{Acknowledgment}
%The author would like to thank\ldots

\appendix
The prior distribution and the measurement likelihood are
\begin{subequations}
\begin{align}
p\left(\xi\right) = & \Npdfbig{\sx}{m}{P} \nonumber \\
& \times \prod_{i=1}^{N_s} \Gammapdf{\gamma_{i}}{\alpha_{i}}{\beta_{i}} \IWishpdf{\ext_{i}}{v_{i}}{V_{i}}, \label{eq:GGIW_prior} \\
p\left(\setZ\left|\xi\right.\right) = & \prod_{i=1}^{N_s} n_{i}! \mathcal{PS}\left(n_{i};\ \gamma_{i}\right) \prod_{j=1}^{n_{i}} \Npdfbig{\sz_{ij}}{H_{i}\sx}{\ext_{i}}. \label{eq:subobject_measurement_likelihood}
\end{align}
\end{subequations}
The problem at hand is to derive the posterior distribution $p\left(\xi|\setZ\right)$ and the corresponding likelihood $\mathcal{L}$,
\begin{align}
p\left(\setZ\left|\xi\right.\right)p\left(\xi\right) = \mathcal{L} p\left(\xi|\setZ\right),
\end{align}
where the posterior $p\left(\xi|\setZ\right)$ is of the same functional form as the prior \eqref{eq:GGIW_prior}, \iep
\begin{align}
p\left(\xi|\setZ\right) = & \Npdfbig{\sx}{m^{+}}{P^{+}} \nonumber \\
& \times \prod_{i=1}^{N_s} \Gammapdf{\gamma_{i}}{\alpha_{i}^{+}}{\beta_{i}^{+}} \IWishpdf{\ext_{i}}{v_{i}^{+}}{V_{i}^{+}}, \label{eq:GGIW_posterior}
\end{align}

The product of Gaussian distributions in \eqref{eq:subobject_measurement_likelihood} can be rewritten as follows
\begin{align}
& \prod_{j=1}^{n_i} \Npdfbig{\sz_{ij}}{H_i\sx}{\ext_i} \label{eq:aux1} \\
= & \left(2\pi\right)^{-n_i d/2}|\ext_i|^{-n_i/2} \nonumber \\
& \times \etr \left(-\frac{1}{2} \left(\sum_{j=1}^{n_i}\left(\sz_{ij}-H_i\sx\right)\left(\sz_{ij}-H_i\sx\right)^{\tp}\right)\ext_i^{-1}\right), \nonumber
\end{align}
where $\etr\left(\cdot\right)=\exp\left(\tr\left(\cdot\right)\right)$ is exponential trace. Define the centroid measurements as
\begin{equation}
\bar{\sz}_i \triangleq \frac{1}{n_i}\sum_{j=1}^{n}\sz_{ij}
\end{equation}
and the scatter matrices as 
\begin{equation}
Z_i\triangleq\sum_{j=1}^{n_i} \left(\sz_{ij} - \bar{\sz}_i\right)\left(\sz_{ij} - \bar{\sz}_i\right)^{\tp},
\end{equation}
and rewrite the summation as
\begin{align}
& \sum_{j=1}^{n_i}\left(\sz_{ij}-{H}_i\sx\right)\left(\sz_{ij}-{H}_i\sx\right)^{\tp} \nonumber \\
= & Z_i+ n_i \left(\bar{\sz}_i - {H}_i\sx\right)\left(\bar{\sz}_i - \tilde{H}_i\sx\right)^{\tp}.
\label{eq:aux2}
\end{align}
Inserting \eqref{eq:aux2} into \eqref{eq:aux1} gives
\begin{subequations}
\begin{align}
&\prod_{j=1}^{n_i} \Npdfbig{\sz_{ij}}{H_i\sx}{\ext_i}\\
=& \left(2\pi\right)^{-\frac{n_i d}{2}}|\ext_i|^{-\frac{n_i}{2}}{ \etr \left(-\frac{1}{2} Z_i\ext_i^{-1}\right)} \nonumber \\
& \times \etr \left(-\frac{1}{2} \left(\bar{\sz}_i - {H}_i\sx\right)\left(\bar{\sz}_i - {H}_i\sx\right)^{\tp}\left(\frac{\ext_i}{n_i}\right)^{-1}\right) \\
=& \left(2\pi\right)^{-\frac{(n_i-1)d}{2}}|\ext_i|^{-\frac{n_i-1}{2}}{n_i}^{-\frac{d}{2}} \nonumber \\
& \times { \etr \left(-\frac{1}{2} Z_i\ext_i^{-1}\right)}\Npdfbig{\bar{\sz}_i}{{H}_i\sx}{\frac{\ext_i}{n_i}} \\
=& \mathcal{A}_{i} \Npdfbig{\bar{\sz}_i}{{H}_i\sx}{\frac{\ext_i}{n_i}}.
\end{align}
\label{eq:aux3}
\end{subequations}
The product of Gaussian distributions in \eqref{eq:subobject_measurement_likelihood} is thus rewritten as
\begin{subequations}
\begin{align}
\prod_{j=1}^{n_i} & \Npdfbig{\sz_{ij}}{H_{i}\sx}{\ext_{i}} = \mathcal{A}_{i} \Npdfbig{\bar{\sz}_{i}}{H_{i}\sx}{\frac{\ext_{i}}{n_{i}}},\\
& \mathcal{A}_{i} = \frac{\left|\ext_{i}\right|^{-\frac{n_i-1}{2}}\etr\left(-\frac{1}{2}Z_{i}\ext_{i}^{-1}\right)}{\left(2\pi\right)^{\frac{(n_i-1)d}{2}} n_i^{\frac{d}{2}} }.
\end{align}%
\end{subequations}%
The product of the likelihood and the prior distribution is
\begin{align}
& p\left(\setZ\left|\xi\right.\right)  p\left(\xi\right) \nonumber \\
= & \Npdfbig{\sx}{m}{P} \left(\prod_{i=1}^{N_s} \mathcal{A}_{i} \Npdfbig{\bar{\sz}_{i}}{H_{i}\sx}{\frac{\ext_{i}}{n_{i}}} \IWishpdf{\ext_{i}}{v_{i}}{V_{i}} \right) \nonumber \\
&\times \left( \prod_{i=1}^{N_s} n_{i}! \mathcal{PS}\left(n_{i};\ \gamma_{i}\right) \Gammapdf{\gamma_{i}}{\alpha_{i}}{\beta_{i}} \right).
\end{align}
For the measurement rates we have the following
\begin{subequations}
\begin{align}
& n_i! \mathcal{PS}\left(n_{i};\ \gamma_{i}\right) \Gammapdf{\gamma_{i}}{\alpha_{i}}{\beta_{i}} \nonumber \\
= & \frac{\beta_{i}^{\alpha_{i}}\gamma_{i}^{\alpha_{i}+n_{i}-1}e^{-(\beta_{i}+1)\gamma_{i}}}{\Gamma\left(\alpha_{i}\right)}  \\
= & \Gammapdf{\gamma_{i}}{\alpha_{i}+n_{i}}{\beta_{i}+1} \frac{\Gamma\left(\alpha_{i}+n_{i}\right)\beta_{i}^{\alpha_{i}}}{\Gamma\left(\alpha_{i}\right)\left(\beta_{i}+1\right)^{\alpha_{i}+n_{i}}} \\
= & \mathcal{L}_{i}^{\gamma} \Gammapdf{\gamma_i}{\alpha_{i}^{+}}{\beta_{i}^{+}}, \label{eq:gamma_meas_update}
\end{align}%
\label{eq:updated_measurement_rates}
\end{subequations}%
where
\begin{subequations}
\begin{align}
\alpha_{i}^{+} &=  \alpha_{i} + n_{i}, \\
\beta_{i}^{+} &= \beta_{i} + 1, \\
\mathcal{L}_{i}^{\gamma} &= \frac{\Gamma\left(\alpha_{i}^{+}\right)\beta_{i}^{\alpha_{i}}}{\Gamma\left(\alpha_{i}\right)\left(\beta_{i}^{+}\right)^{\alpha_{i}^{+}}},
\end{align}
\label{eq:corr_meas_rate}%
\end{subequations}%
The likelihood $\mathcal{L}_{i}^{\gamma}$ is proportional to a negative binomial distribution, see \egp \cite{GelmanCSR:2004}. For the kinematic vector and the random matrices we have
\begin{subequations}
\begin{align}
& \Npdfbig{\sx}{m}{P} \prod_{i=1}^{N_s} \mathcal{A}_{i} \Npdfbig{\bar{\sz}_{i}}{H_{i}\sx}{\frac{\ext_{i}}{n_{i}}} \IWishpdf{\ext_{i}}{v_{i}}{V_{i}} \nonumber \\
= & \Npdfbig{\sx}{m}{P} \Npdfbig{\bar{\sz}}{\mathbb{H}\sx}{\mathbb{X}} \prod_{i=1}^{N_s} \mathcal{A}_{i} \IWishpdf{\ext_{i}}{v_{i}}{V_{i}},
\end{align}
where
\begin{align}
& \bar{\sz} =  \begin{bmatrix} \bar{\sz}_{1}^{\tp} & \cdots & \bar{\sz}_{i}^{\tp} & \cdots & \bar{\sz}_{N_s}^{\tp} \end{bmatrix}^{\tp} \\
& \mathbb{H} =  \begin{bmatrix} H_{1}^{\tp} & \cdots & H_{i}^{\tp} & \cdots & H_{N_s}^{\tp} \end{bmatrix}^{\tp} \\
& \mathbb{X} =  \blkdiag\left( \frac{\ext_{1}}{n_{1}} , \ \ldots , \ \frac{\ext_{i}}{n_{i}} , \ \ldots , \ \frac{\ext_{N_s}}{n_{N_{s}}} \right).
\end{align}%
\end{subequations}%
Using the Kalman filter \cite{Kalman:1960} measurement update we get
\begin{subequations}
\begin{align}
& \Npdfbig{\sx}{m}{P} \Npdfbig{\bar{\sz}}{\mathbb{H}\sx}{\mathbb{X}} \prod_{i=1}^{N_s} \mathcal{A}_{i} \IWishpdf{\ext_{i}}{v_{i}}{V_{i}} \\
= & \Npdfbig{\sx}{\tilde{m}^{+}}{\tilde{P}^{+}} \Npdfbig{\bar{\sz}}{\mathbb{H}m}{\mathbb{H}P\mathbb{H}^{\tp}+\mathbb{X}} \nonumber \\
& \times \prod_{i=1}^{N_s} \mathcal{A}_{i} \IWishpdf{\ext_{i}}{v_{i}}{V_{i}} \label{eq:approximation_equation}
\end{align}
where
\begin{align}
\tilde{m}^{+} = & m+\tilde{K}\left(\bar{\sz}-\mathbb{H}m\right), \\
\tilde{P}^{+} = &  P-\tilde{K}\mathbb{H}P, \\
\tilde{K} = & P\mathbb{H}^{\tp}\tilde{S}^{-1}, \\
\tilde{S} = &  \mathbb{H}P\mathbb{H}^{\tp} + \mathbb{X}. \label{eq:kinematic_correction_innovation_covariance}
\end{align}%
\end{subequations}%
At this point we make two approximations.

\begin{approximation}
In $\tilde{S}$ in \eqref{eq:kinematic_correction_innovation_covariance} the random variable $\mathbb{X}$ is approximated by its expected value
\begin{subequations}
\begin{align}
\hat{\mathbb{X}} = & \E\left[ \mathbb{X} \right] = \blkdiag\left( \frac{\hat{\ext}_{1}}{n_{1}} , \ldots , \frac{\hat{\ext}_{i}}{n_{i}} , \ldots , \frac{\hat{\ext}_{N_s}}{n_{N_{s}}} \right),\\
\hat{\ext}_{i} = & \E\left[\ext_{i}\right] = \frac{V_{i}}{v_{i}-2d-2}
\end{align}%
\end{subequations}%
\hfill$\square$
\end{approximation}
\begin{remark}
This approximation is analogous to an approximation made by Feldmann \etal, see \cite[Equations 31 and 33]{FeldmannFK:2011}.\hfill$\square$
\end{remark}
\begin{approximation}
In $\Npdfbig{\bar{\sz}}{\mathbb{H}m}{\mathbb{H}P\mathbb{H}^{\tp}+\mathbb{X}}$ in \eqref{eq:approximation_equation} the matrix $\mathbb{H}P\mathbb{H}^{\tp}$ is approximated by the block-diagonal matrix
\begin{align}
\blkdiag\left( H_{1}PH_{1}^{\tp} , \ \ldots , \ H_{i}PH_{i}^{\tp} , \ \ldots , \ H_{N_s}PH_{N_s}^{\tp} \right).
\end{align}
\hfill$\square$
\end{approximation}
\begin{remark}
This approximation is necessary to obtain a posterior distribution that is of the same functional form as the prior distribution.
\hfill$\square$
\end{remark}%
Under these approximations instead of \eqref{eq:approximation_equation} we have
\begin{align}
& \Npdfbig{\sx}{m^{+}}{P^{+}} \label{eq:updated_kinematic_vector} \\
\times & \prod_{i=1}^{N_s} \mathcal{A}_{i} \Npdfbig{\bar{\sz}_{i}}{H_{i}m}{H_{i}PH_{i}^{\tp}+\frac{\ext_{i}}{n_{i}}} \IWishpdf{\ext_{i}}{v_{i}}{V_{i}} \nonumber
\end{align}
where
\begin{subequations}
\begin{align}
{m}^{+} = &  m+{K}\left(\bar{\sz}-\mathbb{H}m\right), \\
{P}^{+} = &  P-{K}\mathbb{H}P, \\
{K} = &  P\mathbb{H}^{\tp}{S}^{-1}, \\
{S} = &  \mathbb{H}P\mathbb{H}^{\tp} + \hat{\mathbb{X}}.
\end{align}%
\label{eq:corr_kinematics}
\end{subequations}%
For the factors in the product in \eqref{eq:updated_kinematic_vector}, the Gaussian covariances $H_{i}P H_{i}^{\tp}+\frac{\ext_{i}}{n_{i}}$ can be expanded (\egp using Cholesky Factorization) as
\begin{align}
\left(H_{i}PH_{i}^{\tp}+\frac{\ext_{i}}{n_{i}}\right)^{\frac{1}{2}}\ext_{i}^{-\frac{1}{2}}\ext_{i}\ext_{i}^{-\frac{\tp}{2}}\left(H_{i}PH_{i}^{\tp}+\frac{\ext_{i}}{n_{i}}\right)^{\frac{\tp}{2}}. \label{eq:expanded_Gaussian_covariance}
\end{align}
A third approximation is now made.
\begin{approximation}
Equation \eqref{eq:expanded_Gaussian_covariance} is approximated by
\begin{align}
S_{i}^{\frac{1}{2}} \hat{\ext}_{i}^{-\frac{1}{2}} \ext_{i} \hat{\ext}_{i}^{-\frac{\tp}{2}} S_{i}^{\frac{\tp}{2}}, \qquad S_{i} = H_{i}PH_{i}^{\tp}+\frac{\hat{\ext}_{i}}{n_{i}}.
\end{align}
\hfill$\square$
\end{approximation}
\begin{remark}
This approximation is analogous to an approximation made by Feldmann \etal, see \cite[Equations 38 and 39]{FeldmannFK:2011}.\hfill$\square$
\end{remark}
Under this approximation the Gaussians in \eqref{eq:updated_kinematic_vector} can be rewritten as
\begin{subequations}
\begin{align}
& \Npdfbig{\bar{\sz}_{i}}{H_{i}m}{ S_{i}^{\frac{1}{2}} \hat{\ext}_{i}^{-\frac{1}{2}} \ext_{i} \hat{\ext}_{i}^{-\frac{\tp}{2}} S_{i}^{\frac{\tp}{2}} } \\
= & \left(2\pi\right)^{-\frac{d}{2}} \left| S_{i}^{\frac{1}{2}} \hat{\ext}_{i}^{-\frac{1}{2}} \ext_{i} \hat{\ext}_{i}^{-\frac{\tp}{2}} S_{i}^{\frac{\tp}{2}} \right|^{-\frac{1}{2}} \nonumber \\
& \times \etr\left(-\frac{1}{2}\left(\bar{\sz}_{i}-H_{i}m\right)^{\tp} \left(S_{i}^{\frac{1}{2}} \hat{\ext}_{i}^{-\frac{1}{2}} \ext_{i} \hat{\ext}_{i}^{-\frac{\tp}{2}} S_{i}^{\frac{\tp}{2}}\right)^{-1} \left(\bar{\sz}_{i}-H_{i}m\right) \right) \\
= & \left(2\pi\right)^{-\frac{d}{2}} \left| S_{i}^{\frac{1}{2}} \hat{\ext}_{i}^{-\frac{1}{2}} \right|^{-\frac{1}{2}} \left| \ext_{i} \right|^{-\frac{1}{2}} \left| \hat{\ext}_{i}^{-\frac{\tp}{2}} S_{i}^{\frac{\tp}{2}} \right|^{-\frac{1}{2}} \nonumber \\
& \times \etr\left(-\frac{1}{2}\left(\bar{\sz}_{i}-H_{i}m\right)^{\tp} S_{i}^{-\frac{\tp}{2}} \hat{\ext}_{i}^{\frac{\tp}{2}} \ext_{i}^{-1} \hat{\ext}_{i}^{\frac{1}{2}} S_{i}^{\frac{\-1}{2}} \left(\bar{\sz}_{i}-H_{i}m\right) \right) \\
= & \left(2\pi\right)^{-\frac{d}{2}} \left| \hat{\ext}_{i}^{-\frac{\tp}{2}} S_{i} \hat{\ext}_{i}^{-\frac{1}{2}} \right|^{-\frac{1}{2}} \left| \ext_{i} \right|^{-\frac{1}{2}} \nonumber \\
& \times \etr\left(-\frac{1}{2}\hat{\ext}_{i}^{\frac{1}{2}} S_{i}^{\frac{\-1}{2}} \left(\bar{\sz}_{i}-H_{i}m\right)\left(\bar{\sz}_{i}-H_{i}m\right)^{\tp} S_{i}^{-\frac{\tp}{2}} \hat{\ext}_{i}^{\frac{\tp}{2}} \ext_{i}^{-1} \right) \\
= & \left(2\pi\right)^{-\frac{d}{2}} \left| \hat{\ext}_{i}^{-\frac{\tp}{2}} S_{i} \hat{\ext}_{i}^{-\frac{1}{2}} \right|^{-\frac{1}{2}} \left| \ext_{i} \right|^{-\frac{1}{2}} 
\etr\left(-\frac{1}{2}N_{i} \ext_{i}^{-1} \right)
\end{align}
where
\begin{align}
N_{i} & = \hat{\ext}_{i}^{\frac{1}{2}} S_{i}^{-\frac{1}{2}} \left(\bar{\sz}_{i}-H_{i}m\right)\left(\bar{\sz}_{i}-H_{i}m\right)^{\tp} S_{i}^{-\frac{\tp}{2}} \hat{\ext}_{i}^{\frac{\tp}{2}}
\end{align}
\end{subequations}
The factors in \eqref{eq:updated_kinematic_vector} can now be rewritten as
\begin{subequations}
\begin{align}
\mathcal{A}_{i} & \Npdfbig{\bar{\sz}_{i}}{H_{i}m}{ S_{i}^{\frac{1}{2}} \hat{\ext}_{i}^{-\frac{1}{2}} \ext_{i} \hat{\ext}_{i}^{-\frac{\tp}{2}} S_{i}^{\frac{\tp}{2}} } \IWishpdf{\ext_{i}}{v_{i}}{V_{i}} \\
= & \frac{\left|\ext_{i}\right|^{-\frac{n_i-1}{2}}\etr\left(-\frac{1}{2}Z_{i}\ext_{i}^{-1}\right)}{\left(2\pi\right)^{\frac{(n_i-1)d}{2}} n_i^{\frac{d}{2}} } \nonumber \\
& \times \left(2\pi\right)^{-\frac{d}{2}} \left| \hat{\ext}_{i}^{-\frac{\tp}{2}} S_{i} \hat{\ext}_{i}^{-\frac{1}{2}} \right|^{-\frac{1}{2}} \left| \ext_{i} \right|^{-\frac{1}{2}} 
\etr\left(-\frac{1}{2}N_{i} \ext_{i}^{-1} \right) \nonumber \\
& \times \frac{2^{-\frac{v_{i}-d-1}{2}}\left|V_{i}\right|^{\frac{v_{i}-d-1}{2}}}{\Gamma_{d}\left(\frac{v_{i}-d-1}{2}\right)\left|\ext_{i}\right|^{\frac{v_{i}}{2}}}\etr\left(-\frac{1}{2}\ext_{i}^{-1}V_{i}\right) \\
= & n_i^{-\frac{d}{2}}\left(2 \pi\right)^{-\frac{n_i d}{2}} \left| \hat{\ext}_{i}^{-\frac{\tp}{2}} S_{i} \hat{\ext}_{i}^{-\frac{1}{2}} \right|^{-\frac{1}{2}} \frac{2^{-\frac{v_{i}-d-1}{2}}\left|V_{i}\right|^{\frac{v_{i}-d-1}{2}}}{\Gamma_{d}\left(\frac{v_{i}-d-1}{2}\right)\left|\ext_{i}\right|^{\frac{v_{i}+n_{i}}{2}}} \nonumber \\
& \times  \etr\left(-\frac{1}{2}\left(V_i+Z_i+N_i\right)\ext_{i}^{-1}\right) \\
= & \frac{ n_i^{-\frac{d}{2}}\left(2 \pi\right)^{-\frac{n_i d}{2}}}{ \left| \hat{\ext}_{i}^{-\frac{\tp}{2}} S_{i} \hat{\ext}_{i}^{-\frac{1}{2}} \right|^{\frac{1}{2}}}  \frac{2^{-\frac{v_{i}-d-1}{2}}}{2^{-\frac{v_{i}+n_{i}-d-1}{2}}} \nonumber \\
& \times \frac{\Gamma_{d}\left(\frac{v_{i}+n_i-d-1}{2}\right)}{\Gamma_{d}\left(\frac{v_{i}-d-1}{2}\right)}  \frac{\left|V_{i}\right|^{\frac{v_{i}-d-1}{2}}}{\left|V_{i}+Z_{i}+N_{i}\right|^{\frac{v_{i}+n_i-d-1}{2}}} \nonumber \\
& \times \frac{2^{-\frac{v_{i}+n_{i}-d-1}{2}}\left|V_{i}+Z_{i}+N_{i}\right|^{\frac{v_{i}+n_i-d-1}{2}}}{\Gamma_{d}\left(\frac{v_{i}+n_i-d-1}{2}\right)\left|\ext_{i}\right|^{\frac{v_{i}+n_{i}}{2}}} \nonumber \\
& \times \etr\left(-\frac{1}{2}\left(V_i+Z_i+N_i\right)\ext_{i}^{-1}\right) \\
= & \mathcal{L}^{\sx,\ext}_{i} \IWishpdf{\ext_{i}}{v_{i}^{+}}{V_{i}^{+}}
\end{align}
\label{eq:updated_random_matrices}%
\end{subequations}%
where
\begin{subequations}
\begin{align}
v_{i}^{+} & = v_{i} + n_{i}, \\
V_{i}^{+} & = V_{i} + Z_{i} + N_{i} \\
\mathcal{L}^{\sx,\ext}_{i} & = \frac{ \left(n_i\pi^{n_i}\right)^{-\frac{d}{2}}2^{-\frac{n_i (d-1)}{2}}}{ \left| \hat{\ext}_{i}^{-\frac{\tp}{2}} S_{i} \hat{\ext}_{i}^{-\frac{1}{2}} \right|^{\frac{1}{2}}} \frac{\Gamma_{d}\left(\frac{v_{i}^{+}-d-1}{2}\right)}{\Gamma_{d}\left(\frac{v_{i}-d-1}{2}\right)}  \frac{\left|V_{i}\right|^{\frac{v_{i}-d-1}{2}}}{\left|V_{i}^{+}\right|^{\frac{v_{i}^{+}-d-1}{2}}} \label{eq:kin_ext_likelihood}
\end{align}%
\label{eq:corr_iWish}
\end{subequations}%
The likelihood $\mathcal{L}^{\sx,\ext}_{i}$ is proportional to a generalized matrix variate beta type two distribution, see \egp \cite{GuptaN:2000}.

By the combination of \eqref{eq:updated_measurement_rates}, \eqref{eq:updated_kinematic_vector} and \eqref{eq:updated_random_matrices}, under the three approximations given above, the parameters of the posterior \eqref{eq:GGIW_posterior} are given by \eqref{eq:corr_meas_rate}, \eqref{eq:corr_kinematics}, and \eqref{eq:corr_iWish}, and the likelihood is $\mathcal{L} = \prod_{i=1}^{N_s} \mathcal{L}^{\gamma}_{i} \mathcal{L}^{\sx,\ext}_{i}$, where $\mathcal{L}^{\gamma}_{i}$ is given in \eqref{eq:gamma_meas_update} and $\mathcal{L}^{\sx,\ext}_{i}$ is given in \eqref{eq:kin_ext_likelihood}.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\bibliographystyle{IEEEtran}
\bibliography{PHD_ext_targ_track}

\end{document}

MultipleEllipse.bbl
% Generated by IEEEtran.bst, version: 1.12 (2007/01/11)
\begin{thebibliography}{10}
\providecommand{\url}[1]{#1}
\csname url@samestyle\endcsname
\providecommand{\newblock}{\relax}
\providecommand{\bibinfo}[2]{#2}
\providecommand{\BIBentrySTDinterwordspacing}{\spaceskip=0pt\relax}
\providecommand{\BIBentryALTinterwordstretchfactor}{4}
\providecommand{\BIBentryALTinterwordspacing}{\spaceskip=\fontdimen2\font plus
\BIBentryALTinterwordstretchfactor\fontdimen3\font minus
  \fontdimen4\font\relax}
\providecommand{\BIBforeignlanguage}[2]{{%
\expandafter\ifx\csname l@#1\endcsname\relax
\typeout{** WARNING: IEEEtran.bst: No hyphenation pattern has been}%
\typeout{** loaded for the language `#1'. Using the pattern for}%
\typeout{** the default language instead.}%
\else
\language=\csname l@#1\endcsname
\fi
#2}}
\providecommand{\BIBdecl}{\relax}
\BIBdecl

\bibitem{GilholmGMS:2005}
K.~Gilholm, S.~Godsill, S.~Maskell, and D.~Salmond, ``Poisson models for
  extended target and group tracking,'' in \emph{Proceedings of Signal and Data
  Processing of Small Targets}, vol. 5913.\hskip 1em plus 0.5em minus
  0.4em\relax San Diego, CA, USA: SPIE, Aug. 2005, pp. 230--241.

\bibitem{GilholmS:2005}
K.~Gilholm and D.~Salmond, ``Spatial distribution model for tracking extended
  objects,'' \emph{IEE Proceedings of Radar, Sonar and Navigation}, vol. 152,
  no.~5, pp. 364--371, Oct. 2005.

\bibitem{BaumFH:2012}
M.~Baum, F.~Faion, and U.~D. Hanebeck, ``{Modeling the Target Extent with
  Multiplicative Noise},'' in \emph{Proceedings of the International Conference
  on Information Fusion}, Singapore, Jul. 2012, pp. 2406--2412.

\bibitem{Boersetal:06a}
Y.~Boers, H.~Driessen, J.~Torstensson, M.~Trieb, R.~Karlsson, and
  F.~Gustafsson, ``A track before detect algorithm for tracking extended
  targets,'' \emph{IEE Proceedings Radar, Sonar and Navigation}, vol. 153,
  no.~4, pp. 345--351, Aug. 2006.

\bibitem{PetrovMGA:2011}
N.~Petrov, L.~Mihaylova, A.~Gning, and D.~Angelova, ``A novel sequential monte
  carlo approach for extended object tracking based on border
  parametrisation,'' in \emph{Proceedings of the International Conference on
  Information Fusion}, Chicago, IL, USA, Jul. 2011, pp. 306--313.

\bibitem{Koch:2008}
J.~W. Koch, ``{B}ayesian approach to extended object and cluster tracking using
  random matrices,'' \emph{IEEE Transactions on Aerospace and Electronic
  Systems}, vol.~44, no.~3, pp. 1042--1059, Jul. 2008.

\bibitem{BaumNH:2010}
M.~Baum, B.~Noack, and U.~D. Hanebeck, ``{Extended Object and Group Tracking
  with Elliptic Random Hypersurface Models},'' in \emph{Proceedings of the
  International Conference on Information Fusion}, Edinburgh, UK, Jul. 2010.

\bibitem{GranstromLO:2011}
K.~Gran\-str{\"{o}}m, C.~Lundquist, and U.~Orguner, ``{Tracking Rectangular and
  Elliptical Extended Targets Using Laser Measurements},'' in \emph{Proceedings
  of the International Conference on Information Fusion}, Chicago, IL, USA,
  Jul. 2011, pp. 592--599.

\bibitem{ZhuHL:2011}
H.~Zhu, C.~Han, and C.~Li, ``{An extended target tracking method with random
  finite set observations},'' in \emph{Proceedings of the International
  Conference on Information Fusion}, Chicago, IL, USA, Jul. 2011, pp. 73--78.

\bibitem{ReuterD:2011}
S.~Reuter and K.~Dietmayer, ``Pedestrian tracking using random finite sets,''
  in \emph{Proceedings of the International Conference on Information Fusion},
  Chicago, IL, USA, Jul. 2011, pp. 1101--1108.

\bibitem{DegermanWS:2011}
J.~Degerman, J.~Wintenby, and D.~Svensson, ``Extended target tracking using
  principal components,'' in \emph{Proceedings of the International Conference
  on Information Fusion}, Chicago, IL, USA, Jul. 2011, pp. 330--337.

\bibitem{LanRL:2012}
J.~Lan and X.-R. Li, ``Tracking of extended object or target group using random
  matrix -- part {I}: New model and approach,'' in \emph{Proceedings of the
  International Conference on Information Fusion}, Singapore, Jul. 2012, pp.
  2177--2184.

\bibitem{ReuterWD:2012}
S.~Reuter, B.~Wilking, and K.~Dietmayer, ``Methods to model the motion of
  extended objects in multi-object {B}ayes filters,'' in \emph{Proceedings of
  the International Conference on Information Fusion}, Singapore, Jul. 2012,
  pp. 527--534.

\bibitem{LundquistGO:11}
C.~Lundquist, K.~Gran\-str{\"{o}}m, and U.~Orguner, ``{Estimating the Shape of
  Targets with a PHD Filter},'' in \emph{Proceedings of the International
  Conference on Information Fusion}, Chicago, IL, USA, Jul. 2011, pp. 49--56.

\bibitem{BaumH:2011}
M.~Baum and U.~D. Hanebeck, ``{Shape Tracking of Extended Objects and Group
  Targets with Star-Convex RHMs},'' in \emph{Proceedings of the International
  Conference on Information Fusion}, Chicago, IL, USA, Jul. 2011, pp. 338--345.

\bibitem{LanRL:2014}
J.~Lan and X.~Rong-Li, ``Tracking maneuvering non-ellipsoidal extended object
  or target group using random matrix,'' \emph{IEEE Transactions on Signal
  Processing}, 2014.

\bibitem{LanRL:2012irregular}
J.~Lan and X.-R. Li, ``Tracking of extended object or target group using random
  matrix -- part {II}: Irregular object,'' in \emph{Proceedings of the
  International Conference on Information Fusion}, Singapore, Jul. 2012, pp.
  2185--2192.

\bibitem{mahler_FUSION_2009_extTarg}
R.~Mahler, ``{PHD} filters for nonstandard targets, {I:} {E}xtended targets,''
  in \emph{Proceedings of the International Conference on Information Fusion},
  Seattle, WA, USA, Jul. 2009, pp. 915--921.

\bibitem{GranstromLO:2010}
K.~Gran\-str{\"{o}}m, C.~Lundquist, and U.~Orguner, ``A {G}aussian mixture
  {PHD} filter for extended target tracking,'' in \emph{Proceedings of the
  International Conference on Information Fusion}, Edinburgh, UK, Jul. 2010.

\bibitem{GranstromLO:2012}
------, ``{Extended Target Tracking using a Gaussian Mixture PHD filter},''
  \emph{IEEE Transactions on Aerospace and Electronic Systems}, vol.~48, no.~4,
  pp. 3268--3286, Oct. 2012.

\bibitem{GranstromO:2012a}
K.~Gran\-str{\"{o}}m and U.~Orguner, ``{A PHD filter for tracking multiple
  extended targets using random matrices},'' \emph{IEEE Transactions on Signal
  Processing}, vol.~60, no.~11, pp. 5657--5671, Nov. 2012.

\bibitem{LundquistGO:2013}
C.~Lundquist, K.~Gran\-str{\"{o}}m, and U.~Orguner, ``An extended target {CPHD}
  filter and a gamma {G}aussian inverse {W}ishart implementation,'' \emph{IEEE
  Journal of Selected Topics in Signal Processing, Special Issue on
  Multi-target Tracking}, vol.~7, no.~3, pp. 472--483, Jun. 2013.

\bibitem{SwainC:2010}
A.~Swain and D.~Clark, ``Extended object filtering using spatial independent
  cluster processes,'' in \emph{Proceedings of the International Conference on
  Information Fusion}, Edinburgh, UK, Jul. 2010.

\bibitem{SwainC:2012}
------, ``The {PHD} filter for extended target tracking with estimable shape
  parameters of varying size,'' in \emph{Proceedings of the International
  Conference on Information Fusion}, Singapore, Jul. 2012.

\bibitem{GuptaN:2000}
A.~K. Gupta and D.~K. Nagar, \emph{Matrix variate distributions}, ser. Chapman
  \& Hall/CRC monographs and surveys in pure and applied mathematics.\hskip 1em
  plus 0.5em minus 0.4em\relax Chapman \& Hall, 2000.

\bibitem{GranstromO:2014}
K.~Gran\-str{\"{o}}m and U.~Orguner, ``{A New Prediction Update for Extended
  Target Tracking with Random Matrices},'' \emph{IEEE Transactions on Aerospace
  and Electronic Systems}, vol.~50, no.~2, Apr. 2014.

\bibitem{FeldmannF:2008}
M.~Feldmann and D.~Fr{\"{a}}nken, ``{Tracking of Extended Objects and Group
  Targets using Random Matrices - A New Approach},'' in \emph{Proceedings of
  the International Conference on Information Fusion}, Cologne, Germany, Jul.
  2008.

\bibitem{FeldmannFK:2011}
M.~Feldmann, D.~Fr{\"{a}}nken, and J.~W. Koch, ``Tracking of extended objects
  and group targets using random matrices,'' \emph{IEEE Transactions on Signal
  Processing}, vol.~59, no.~4, pp. 1409--1420, Apr. 2011.

\bibitem{Orguner:2012}
U.~Orguner, ``A variational measurement update for extended target tracking
  with random matrices,'' \emph{IEEE Transactions on Signal Processing},
  vol.~60, no.~7, pp. 3827--3834, Jul. 2012.

\bibitem{LianHLYZ:2010}
F.~Lian, C.-Z. Han, W.-F. Liu, X.-X. Yan, and H.-Y. Zhou, ``Sequential {M}onte
  {C}arlo implementation and state extraction of the group probability
  hypothesis density filter for partly unresolvable group targets-tracking
  problem,'' \emph{IET Radar, Sonar and Navigation}, vol.~4, no.~5, pp.
  685--702, Oct. 2010.

\bibitem{StreitL:1993}
R.~Streit and T.~E. Luginbuhl, ``A probabilistic multihypothesis tracking
  algorithm without enumeration and pruning,'' in \emph{In proceedings of the
  Sixth Joint Service Data Fusion Symposium}, Laurel, MD, USA, Jun. 1993, pp.
  1015--1024.

\bibitem{WienekeK:2010}
W.~Wieneke and J.~W. Koch, ``Probabilistic tracking of multiple extended
  targets using random matrices,'' in \emph{Proceedings of SPIE Signal and Data
  Processing of Small Targets}, Orlando, FL, USA, Apr. 2010.

\bibitem{WienekeD:2011}
M.~Wieneke and S.~Davey, ``Histogram pmht with target extent estimates based on
  random matrices,'' in \emph{Proceedings of the International Conference on
  Information Fusion}, Chicago, IL, USA, Jul. 2011, pp. 1--8.

\bibitem{WienekeK:2012}
M.~Wieneke and W.~Koch, ``A {PMHT} approach for extended objects and object
  groups,'' \emph{IEEE Transactions on Aerospace and Electronic Systems},
  vol.~48, no.~3, pp. 2349--2370, 2012.

\bibitem{GranstromO:2012c}
K.~Gran\-str{\"{o}}m and U.~Orguner, ``{Estimation and Maintenance of
  Measurement Rates for Multiple Extended Target Tracking},'' in
  \emph{Proceedings of the International Conference on Information Fusion},
  Singapore, Jul. 2012, pp. 2170--2176.

\bibitem{BarShalomWT:2011}
Y.~Bar-Shalom, P.~K. Willett, and X.~Tian, \emph{Tracking and data fusion, a
  handbook of algorithms}.\hskip 1em plus 0.5em minus 0.4em\relax YBS, 2011.

\bibitem{RongLiJ:2003}
X.-R. Li and V.~Jilkov, ``Survey of maneuvering target tracking: {P}art {I}.
  {D}ynamic models,'' \emph{IEEE Transactions on Aerospace and Electronic
  Systems}, vol.~39, no.~4, pp. 1333--1364, Oct. 2003.

\bibitem{Kalman:1960}
R.~E. Kalman, ``{A New Approach to Linear Filtering and Prediction Problems},''
  \emph{Transactions of the ASME - Journal of Basic Engineering}, vol.~82, no.
  Series D, pp. 35--45, Mar. 1960.

\bibitem{Jazwinski:1970}
A.~H. Jazwinski, \emph{{Stochastic Processes and Filtering Theory}}.\hskip 1em
  plus 0.5em minus 0.4em\relax Academic Press, 1970.

\bibitem{Bishop:2006}
C.~M. Bishop, \emph{Pattern recognition and machine learning}.\hskip 1em plus
  0.5em minus 0.4em\relax New York, USA: Springer, 2006.

\bibitem{HastieTF:2009}
T.~Hastie, R.~Tibshirani, and J.~H. Friedman, \emph{The elements of statistical
  learning : {D}ata mining, inference, and prediction}, 2nd~ed.\hskip 1em plus
  0.5em minus 0.4em\relax New York, NY, USA: Springer, 2009.

\bibitem{DempsterLR:1977}
A.~P. Dempster, N.~M. Laird, and D.~B. Rubin, ``Maximum likelihood from
  incomplete data via the {EM} algorithm,'' \emph{Journal of the Royal
  Statistical Society, Series B (Methodological)}, vol.~39, no.~1, pp. 1--38,
  1977.

\bibitem{GranstromO:2012e}
K.~Gran\-str{\"{o}}m and U.~Orguner, ``{On the Reduction of Gaussian inverse
  Wishart mixtures},'' in \emph{Proceedings of the International Conference on
  Information Fusion}, Singapore, Jul. 2012, pp. 2162--2169.

\bibitem{GelmanCSR:2004}
A.~Gelman, J.~B. Carlin, H.~S. Stern, and D.~B. Rubin, \emph{Bayesian Data
  Analysis}, ser. Texts in Statistical Science.\hskip 1em plus 0.5em minus
  0.4em\relax Chapman \& Hall/CRC, 2004.

\end{thebibliography}
