Least square ellipsoid fitting using iterative orthogonal transformations

TL;DR

This paper introduces a numerically stable iterative method for fitting ellipsoids to minimal data points, accurately determining shape, orientation, and axes, with applications in gravitational wave data analysis.

Contribution

The paper presents a novel, stable ellipsoid fitting algorithm that works with minimal data points and accurately retrieves shape and orientation parameters.

Findings

Effective on simulated data sets

Applicable to highly elongated and arbitrarily oriented ellipsoids

Potential use in gravitational wave data analysis

Abstract

We describe a generalised method for ellipsoid fitting against a minimum set of data points. The proposed method is numerically stable and applies to a wide range of ellipsoidal shapes, including highly elongated and arbitrarily oriented ellipsoids. This new method also provides for the retrieval of rotational angle and length of semi-axes of the fitted ellipsoids accurately. We demonstrate the efficacy of this algorithm on simulated data sets and also indicate its potential use in gravitational wave data analysis.

Figures

Tables

Table 1. TABLE I: Data Set 1 ( χ = 1.5 𝜒 1.5 \chi=1.5 , aligned)

Input	This work		Li(2004)	Ying(2012)
	Fisher matrix	Random matrix
$A=12.0$	$12.0$	$12.0$	$12.0$	$12.0$
$B=10.0$	$10.0$	$10.0$	$10.0$	$10.0$
$C=8.0$	$8.0$	$8.0$	$8.0$	$8.0$
$\alpha =0.0,\beta =0.0,\gamma =0.0$	$\alpha =0.0,\beta =0.0,\gamma =0.0$	$\alpha =0.0,\beta =0.0,\gamma =0.0$	-	-

Table 2. TABLE II: Data Set 1 ( χ = 5 𝜒 5 \chi=5 , aligned)

Input	This work		Li(2004)	Ying(2012)
	Fisher matrix	Random matrix
$A=5.0$	$5.0$	$5.0$	$4.086188652$	$5.0$
$B=3.0$	$3.0$	$3.0$	$3.586034787$	$3.0$
$C=1.0$	$1.0$	$1.0$	$2.684837382$	$1.0$
$\alpha =0.0,\beta =0.0,\gamma =0.0$	$\alpha =0.0,\beta =0.0,\gamma =0.0$	$\alpha =0.0,\beta =0.0,\gamma =0.0$	-	-

Table 3. TABLE III: Data Set 3 ( χ = 10 𝜒 10 \chi=10 , aligned)

Input	This work		Li(2004)	Ying(2012)
	Fischer matrix	Random matrix
$A=10.0$	$10.0$	$10.0$	$10.050239317$	$10.0$
$B=6.0$	$6.0$	$6.0$	$5.103102823$	$6.0$
$C=1.0$	$1.0$	$1.0$	$4.547883756$	$1.0$
$\alpha =0.0,\beta =0.0,\gamma =0.0$	$\alpha =0.0,\beta =0.0,\gamma =0.0$	$\alpha =0.0,\beta =0.0,\gamma =0.0$	-	-

Table 4. TABLE IV: Data Set 4 ( χ = 1.5 𝜒 1.5 \chi=1.5 , non-aligned)

Input	This work		Li(2004)	Ying (2012)
	Fisher matrix	Random matrix
$A=12.0$	$12.0$	$12.0$	$12.0$	$16.4749$
$B=10.0$	$10.0$	$10.0$	$10.0$	$9.2354$
$C=8.0$	$8.0$	$8.0$	$8.0$	$7.6336$
$\alpha =30$	$30.0$	$30.0$	-	-
$\beta =80$	$80.0$	$80.0$	-	-
$\gamma =70$	$70.0$	$70.0$	-	-

Table 5. TABLE V: Data Set 5 ( χ = 5 𝜒 5 \chi=5 , non-aligned)

Input	This work		Li(2004)	Ying(2012)
	Fisher matrix	Random matrix
$A=1.0$	$0.999257113$	$0.999128747$	$1.280262088$	$0.9129$
$B=3.0$	$3.002419519$	$3.002838810$	$2.550536762$	$1.2791$
$C=5.0$	$5.000012418$	$5.000014725$	$5.058022788$	$1.3335$
$\alpha =70$	$69.980565$	$69.977193$	-	-
$\beta =10$	$9.989666$	$9.987882$	-	-
$\gamma =30$	$30.007390$	$30.008666$	-	-

Table 6. TABLE VI: Data Set 6 ( χ = 10 𝜒 10 \chi=10 , non-aligned)

Input	Our Method		Li(2004)	Ying(2012)
	Fisher matrix	Random matrix
$A=10.0$	$10.0$	$10.0$	$5.076706206$	$1.9486$
$B=3.0$	$3.0$	$3.0$	$2.862531749$	$1.5969$
$C=1.0$	$1.0$	$1.0$	$2.492054846$	$0.8945$
$\alpha =50$	$50.0$	$50.0$	-	-
$\beta =60$	$60.0$	$60.0$	-	-
$\gamma =40$	$40.0$	$40.0$	-	-

Table 7. TABLE VII: Reconstructed semi-axes and Euler angles for constant match ellipsoid as described in Section V

$A=2.14$

$B=0.047$

$C=0.004$

$\alpha =33.72$

$\beta =7.72$

$\gamma =19.53$

Equations

$$a{x}^2+b{y}^2+c{z}^2+2fyz+2gxz+2hxy+2px+2qy+2rz+d=0.$$

$$kJ-I^2=1.$$

$$I\equiv a+b+c,$$

$$J\equiv ab+bc+ac-f^2-g^2-h^2,$$

$${\mathbf{X}}_{\mathbf{i}}=(x_i^2,y_i^2,z_i^2,2y_i{z}_i,2x_i{z}_i,2x_i{y}_i,2x_i,2y_i,2z_i,1{)}^T.$$

$$\Omega =\sum _{i=1}^N({\mathbf{v}}^T{\mathbf{X}}_{\mathbf{i}}{)}^2$$

$$\mathbf{v}\equiv (a,b,c,f,g,h,p,q,r,d{)}^T$$

$${\mathbf{D}}^T\mathbf{v}=0,$$

$$\mathbf{C}=\left[\begin{array}{cccccccccc}-1& \frac{k}{2}-1& \frac{k}{2}-1& 0& 0& 0& 0& 0& 0& 0\\ \frac{k}{2}-1& -1& \frac{k}{2}-1& 0& 0& 0& 0& 0& 0& 0\\ \frac{k}{2}-1& \frac{k}{2}-1& -1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& -k& 0& 0& 0& 0& 0& 0\\ .0& 0& 0& 0& -k& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& -k& 0& 0& 0& 0\end{array}\right].$$

$$\mathcal{L}(\mathbf{v},\lambda )=\Omega -\lambda ({\mathbf{v}}^T\mathbf{Cv}-1),$$

$${\mathbf{DD}}^T\mathbf{v}=\lambda \mathbf{Cv},$$

$${\mathbf{v}}^T\mathbf{Cv}=1,$$

$${\mathbf{A}}^{\mathbf{T}}\mathbf{KA}=-d$$

$$\mathbf{K}=\left[\begin{array}{cccc}a& h& g& \\ h& b& f& \\ g& f& c& \end{array}\right],$$

$$\mathbf{A}=\left[\begin{array}{ccc}x& y& z\end{array}\right].$$

$$\mathbf{R}=\left[\begin{array}{ccc}{R}_{11}& R_{12}& R_{13}\\ R_{21}& R_{22}& R_{23}\\ R_{31}& R_{32}& R_{33}\end{array}\right]$$

$$\mathbf{R}={\mathbf{R}}_{\mathbf{Z}}(\gamma ){\mathbf{R}}_{\mathbf{Y}}(\beta ){\mathbf{R}}_{\mathbf{X}}(\alpha ),$$

$$\{p_i\}={\mathbf{R}}^T\left[\begin{array}{ccc}A\cos {\theta }_i\cos {\varphi }_i& B\cos {\theta }_i\sin {\varphi }_i& C\sin {\theta }_i\end{array}\right]$$

$$\begin{array}{rl}{R}_{11}& =\cos \alpha \cos \beta ,\\ R_{12}& =\sin \gamma \sin \beta \cos \alpha -\cos \gamma \sin \alpha ,\\ R_{13}& =\cos \gamma \sin \beta \cos \alpha +\sin \gamma \sin \alpha ,\\ R_{21}& =\cos \beta \sin \alpha ,\\ R_{22}& =\sin \gamma \sin \alpha \sin \beta +\cos \gamma \cos \alpha ,\\ R_{23}& =\cos \gamma \sin \alpha \sin \beta -\sin \gamma \cos \alpha ,\\ R_{31}& =-\sin \beta ,\\ R_{32}& =\sin \gamma \cos \beta ,\\ R_{33}& =\cos \gamma \cos \beta .\end{array}$$

$$⟨h({\vec{\lambda }}_1),h({\vec{\lambda }}_2)⟩\equiv 4\mathrm{Re}{\int }_0^{\infty }df\frac{{\tilde{h}}^{\ast }(f;{\vec{\lambda }}_1)\tilde{h}(f;{\vec{\lambda }}_2)}{S_h(f)},$$

$$M(\vec{\lambda},\Delta{\vec{\lambda}})\approx 1+\sum_{i,j=1}^{N}\frac{1}{2}\left(\frac{\partial^{2}M}{\partial\Delta{\lambda^{i}}\partial\Delta{\lambda^{j}}}\right)\biggr{\rvert}_{\Delta\vec{\lambda}=0}\Delta{\lambda}^{i}\Delta{\lambda}^{j}.$$

$$\Delta s^2$$

$$=$$

$$g_{ij}=-\frac{1}{2}\left(\frac{\partial^{2}M}{\partial\Delta{\lambda^{i}}\Delta{\lambda}^{j}}\right)\biggr{\rvert}_{\Delta{\vec{\lambda}}=0}.$$

$$g_{ab}=⟨\partial h/\partial {\lambda }_a,\partial h/\partial {\lambda }_b⟩;a,b=1,2,3;{\lambda }_{a,b}\in \vec{\lambda }.$$

Full text

I Introduction

Reconstructing 3-dimensional (3D) ellipsoidal surfaces from discrete data points is a well-studied problem in the field of computer vision, pattern recognition, astronomy and medical image processing. Such 3D ellipsoid models find use in a broad class of applications: (a) One of the primary goals of computer vision is to find a suitable 3D shape descriptor for arbitrary shapes while preserving geometrical information as much as possible. Existing shape representations based on spherical descriptors are limited in scope as they are appropriate only for nearly spherical bodies. On the other hand, the ellipsoidal descriptors provide a closer approximation to irregular 3D shapes. Khatun et al. [1] have proposed a method where ellipsoidal 3D shape representations have been used to retrieve arbitrary 3D shapes. (b) Gait analysis is utilised for a systematic study of animal locomotion by measuring and tracking their body mechanics. Gathering gait features and interpretation of the gait dynamics is a challenging task as current techniques are computationally expensive. 3D ellipsoid fitting methods have been found useful in this line of research. Sivapalan et al. [2] have proposed a fast 3D ellipsoid based gait recognition algorithm using a 3D voxel model, (c) galaxies are often modelled as 3-dimensional ellipsoids whose parameters are determined from images recorded by telescopes. Compère et al. [3] have recently proposed a three-dimensional galaxy fitting algorithm to extract parameters of the bulge, long bar, disc and a central point source from broadband images of galaxies.

Several techniques exist for fitting ellipsoids to a set of data points and can be broadly classified into projection based algorithms [4, 5] and nonlinear optimisation based surface fitting algorithms [6, 7]. Projection based fitting algorithms are also further organised into two categories: namely, orthographic and line integral based.

The basic idea in orthographic projection is to use matrix operators to project a 3D shape onto planes. For the case of a 3D ellipsoid, three different orthographic transformations are possible along three orthogonal planes, and projected shapes are 2D ellipses. If the parameters of the projected ellipses are deciphered, then 3D rotation between two successive projections can be detected after characterising the variations of the semi-axes length and the orientation of the projected ellipses. On the other hand, the line integral projection based methods commonly use projection contours to reconstruct ellipsoids. The general second-degree equation of an ellipsoid is used to construct the line integral projection model. One of the shortcomings of such methods is that one needs prior information about the projected ellipses to determine the angle and axis of rotation accurately.

In nonlinear optimisation techniques for fitting ellipsoidal surfaces, the latter is modelled as a bounded surface through a family of polynomials which are then fitted using standard nonlinear optimisation methods. The problem with such techniques is that due to high non-linearity of the model, the optimal solution may get stuck in local solutions leaving the resulting surface unbounded. Therefore the proper solution can not guarantee closed bounded solution of the desired surface. To overcome this problem Li et al. [8] have prescribed an algorithm to obtain closed form of the resulting ellipsoidal surface by providing an additional constraint. The fitting algorithm works robustly for ellipsoids whose short radius is at least half of their major radius and for which the semi-axes of the model ellipsoid are aligned along the co-ordinates. More recently, Ying et al. (2012) [9] have proposed a least-square ellipsoid fitting algorithm by extending the 2-dimensional ellipse fitting algorithm given by Fitzgibbon et al. [10].

In this work, we propose a stable algorithm that can fit an ellipsoidal surface to a given set of data points and can detect the rotational angle as well as semi-axes length with significant improvement over Li (2004) [8]. The new method is applicable even to extreme cases where the ellipsoid is highly elongated and arbitrarily oriented to a rigid frame of reference. Our primary motivation is to extend their idea in such an algorithmic form, which can produce the best fit for any kind of the ellipsoidal surface. Also, we describe a method for the retrieval of the orientation of such ellipsoids without assuming any prior information.

This paper is organised as follows: In Section-II, we present a concise description of Li et al. [8], establishing the notation used in this article and highlighting the salient features of their algorithm. Section-III describes our proposed method based on the general equation of an ellipsoid. In section-IV, we describe the algorithm for retrieval of the orientation of the reconstructed ellipsoid, followed by a demonstration of the efficacy of our method using synthetic data. We then present a case study where this approach is applied to the field of gravitational wave data analysis in Section-V.

Finally, we make some general comments on the results obtained in this paper.

II Previous Work:

The general equation of the second degree in three variables $(x,y,z)$ representing a conic is given by:

[TABLE]

As shown in [8], Eq. (1) represents an ellipsoid under the constraint

[TABLE]

where

[TABLE]

[TABLE]

and $k$ is a positive number. For ellipsoids with comparable semi-axes lengths, $k\sim 4$.

Let P = $\mathbf{p_{i}(x_{i},y_{i},z_{i})},\{$i = 1,2,..,N$\}$ be the coordinates of $N$ points with respect to a fixed frame of reference $XYZ$ (refer Figure 1b) to which an ellipsoid is to be fitted. Further, let the ellipsoid be arbitrarily oriented in this frame. For every point $\mathbf{p}_{i}$, one defines a column $\mathbf{X_{i}}$ of the design matrix $\mathbf{D}$ as:

[TABLE]

To fit an ellipsoidal surface, each data point must satisfy the quadratic Eq. (1) with the constraint defined in Eq. (2). The algebraic distance $\Omega$ between the model and the set of data points defined as,

[TABLE]

must be minimized with respect to $\mathbf{v}$ in order to find the best fit, where

[TABLE]

is the set of unknown parameters whose values are to be determined. Therefore the ellipsoid fitting problem can be mapped to an optimisation problem that can be solved using standard least square methods.

It is obvious that Eq. (1) can be written in matrix form as the following system of linear equations:

[TABLE]

in terms of the design matrix $\mathbf{D}=(X_{1},X_{2},...,X_{i})$ of order $10\times N$, where $N\geq 10$. The geometric distance above can also be written in matrix form as $\Omega=||{\mathbf{Dv}}||^{2}=\mathbf{v}^{T}\mathbf{D}^{T}\mathbf{vD}$, which is to be minimized subject to the constraint given in Eq. (2). The latter can also be written in matrix form as $\mathbf{v}^{T}\mathbf{Cv}=1$, where

[TABLE]

The Lagrangian of this optimization problem is defined as,

[TABLE]

where $\lambda$ is the scalar Lagrange multiplier. Using the standard Lagrange multiplier method [11], we set ${\partial\mathcal{L}}/{\partial\bf v}=0$ and ${\partial\mathcal{L}}/{\partial\lambda}=0$ leading to

[TABLE]

and

[TABLE]

respectively.

Eq. (11) is in the form of a generalized eigenvalue equation which can be solved for $\lambda$ and $\bf{v}$. The eigenvectors and corresponding eigenvalues can be used to determine the semi-axes length and its orientation as explained later in this paper. The optimal value of the parameter $k$ in Eq. (2) depends on the input data set: for a given data set; the optimisation leads to the correct value of $k$ which is to be determined iteratively.

Ying et al. (2012) [9] have also established a similar type of generalised eigenvalue equation by introducing the concept of a random plane that intersects the quadratic surface defined in Eq. (1). Under this paradigm, a quadratic surface is deemed to be an ellipsoid if its intersection with any random plane is an ellipse. Incorporating this idea leads to some modifications to the constraint matrix for the generalised eigenvalue system in Eq. (11). The eigenvector corresponding to a unique positive eigenvalue provides the required solution.

III Methodology

Eq. (1) can be normalized by $d(\neq 0)$. Therefore the actual number of unknown parameters involved in the system is 9. Further, if the centre of the ellipsoid is known, one can fix it at the origin $(0,0,0)$ without any loss of generality, in which case, the number of unknown parameters further reduces to 6. Thus a minimum of six unique data points is sufficient to find $\mathbf{v}$ unambiguously.

Eq. (1) can be written in matrix form as:

[TABLE]

where

[TABLE]

[TABLE]

The Fisher information matrix $\mathbf{K}$ is constructed from the elements of $\mathbf{v}$, and our aim in fitting the ellipsoid is to reconstruct this matrix robustly from the given data points. The eigenvectors of $\bf{K}$ are aligned along the principal directions of the ellipsoid. Off-diagonal terms signify cross-correlation among the variables and allude to the fact that these axes are not aligned along the rigid frame of reference $XYZ$. Conversely, $f=g=h=0$ implies that principle axes of the ellipsoid are aligned along the rigid frame of reference.

Starting from an initial estimate, we aim to find a conformal transformation through a rotation matrix $\mathbf{R}$ in an iterative fashion, in which $\bf K$ becomes diagonal. Note that rotation matrices must be orthogonal, i.e. $\mathbf{R}\mathbf{R^{T}}=\mathbf{I}$.

An initial estimate of $\mathbf{K}$ can be made by using the fact that the Fisher information matrix is equal to the inverse of the data covariance matrix. Uniform sampling of data points over the ellipsoid can lead to a good initial estimate of $\bf K$ by this method, but pathological cases may arise when all the data points are sampled from a narrow region on the ellipsoidal surface. In the latter case, the inverse of the data covariance matrix (if it exists) may not be a good initial estimate of $\bf{K}$. The corresponding $\mathbf{R}$ is constructed from the eigenvectors of $\mathbf{K}$.

Alternatively, $\mathbf{R}$ can be constructed from a $\it{random}$ positive definite matrix $\mathbf{G}\in\mathcal{N}(0,1)$ in such a way that each column of $\mathbf{R}$ is formed from the eigenvectors of the covariance of $\mathbf{G}$.

Regardless of the method used to initialize $\mathbf{R}$, we show in Fig 2a-2b that our algorithm achieves the convergence criteria leading to the diagonal form of $\bf{K}$ after a few iterative steps. Since no prior information is assumed about the way in which data points are sampled from the surface, it is advisable to initialise $\bf R$ using the random matrix method.

In a single iteration of this algorithm, $\mathbf{R}$ is used to project the data points which are then used to find the best fitting ellipsoid using the method of least squares as outlined in Section II. This process leads to the best fit $\mathbf{K}$ whose eigenvectors are used to further refine $\mathbf{R}$. This process continues until the desired termination criteria are met. In the process, the fixed axes $XYZ$ undergoes a series of successive conformal transformations until it aligns with the principal directions of the ellipsoid. The successive updates to $\mathbf{R}$ are recorded and used to reconstruct $\mathbf{K}$ in $XYZ$ by applying an inverse transformation: $\mathbf{K}_{XYZ}=\mathbf{R^{T}}\mathbf{KR}$.

The steps of this method are given in Algorithm 1 and illustrated in Figure 1.

III-A Convergence criteria

As stated above, the iterative algorithm proceeds by projecting the data points through a succession of conformal transformations until $\mathbf{K}$ becomes diagonal. At convergence, the off-diagonal terms $f,g,h$ become nearly zero, and the trace of the eigenvector matrix111Eigen vector matrix refers to the matrix whose columns are the eigenvectors in a decreasing order of the eigenvalues of $\mathbf{K}$ is nearly equal to $3$.

The initial estimate of $\mathbf{R}$ may be quite off the mark in which case; the least square fitting algorithm may require an enormous value of $k$ in the constraint equation Eq. (2) for meeting the least-square convergence criteria. It can even lead to a failure to find the optimal least square solution. This can be solved by restricting $k$: $4\leq k\leq k_{max}$ to a maximum upper limit $k_{max}\sim 10^{10}$. As a consequence, the least square ellipsoid fittings in the early stages of the iterative algorithm may be sub-optimal, but it does not affect the outcome. As successive conformal transformations are applied, one requires progressively smaller values of $k$ for convergence at every iteration. Restricting $k$ to an upper limit also makes it computationally more efficient.

If the ellipsoid is aligned along fixed frame, then an initial choice of $\bf R=\mathcal{I}_{3\times 3}$ will further reduce the number of iterative transformations to a one step fitting problem.

III-B Efficient retrieval of orientation angles

A general Euler rotation matrix is of the following form.

[TABLE]

This matrix can be factorized as a product of rotations in the following sequence:

[TABLE]

where $\alpha$, $\beta$, $\gamma$ represents the Euler angles corresponding to $X$, $Y$, and $Z$ axes respectively. The elements $R_{ij}$ of $\mathbf{R}$ represent a specific function of the Euler angles as given in Eq. (19). These functions can be inverted to solve for $\alpha$, $\beta$, $\gamma$. Note that these are an overdetermined set of 9 equation in 3 unknowns. We may use any three to solve for the Euler angles and use the others for consistency check.

As we do not assume any prior information about the orientation and shape of the ellipsoid, we initialize $\mathbf{R}$ randomly which undergoes a series of refinements through the iterative fitting procedure. At the end of this fit, the fixed frame $XYZ$ when transformed through $\mathbf{R}$ aligns itself along the ellipsoid axes. However, this alignment can happen in 6 different ways. This is easily understood from Figure 1a - the transformed X-axis can align along $\pm X_{1}$ or $\pm Y_{1}$ or $\pm Z_{1}$ directions. The convention used for determining the Euler angles [12] is such that the transformed $X$ axis lies along the principal $X_{1}$ direction of the ellipsoid. Therefore, as highlighted in Algorithm 2, one must scan through all 6 column permutations of the rotation matrix $\mathbf{R}$ to re-orientate the axes correctly.

III-C Computational complexity

The proposed iterative method depends on two different convergence criteria (lsqConv (Eq. (2)) and conv (III-A)). In the inner loop, the convergence criteria lsqConv optimises the value of $k$. The time complexity of this circuit depends on solving the generalised eigenvalue problem as given by Eq. (11). On the other hand, the conv criteria in the outer loop optimise the alignment of the principal axes of the fitted ellipsoid.

From Algorithm 1, the total time complexity for each iteration in the inner loop can be seen to be $\mathcal{O}(10\;N^{2}+10^{3})$, where $N$ is the number of supplied data points. The number of iterations depends on $\chi$. In our numerical experiments, we have noted that extremely elongated ellipsoids (large $\chi$) can take 35 – 40 iterations (with six input data points). On the other hand for smaller values of $\chi$, only a few iterations are enough.

It is also evident that every iteration of the outer loop requires $\mathcal{O}(3N^{2}+3^{3})$ operations from the matrix multiplications in Step $4$ and $15$ respectively. Just like the inner loop, the number of iterations for the outer loop also depend on $\chi$. In our experience with simulated data sets, this can be as high as 15 – 20 for extremely elongated ellipsoids. For input data from moderately elongated ellipsoids, the outer loop is seen to converge within 2 – 4 iterations.

The proposed algorithm is robust even for a small number of data points: results presented in this paper are for the minimum set of data points.

IV Implementation and Experimental Results

In this section, we describe the generation of synthetic data points to confront it against the new algorithm and present the results. We also compare these results against Li’s (2004) [8] and Ying’s (2012) [9] method wherever possible.

Synthetic Data set

Let $\alpha$, $\beta$, $\gamma$ be the Euler angles and $A$, $B$, $C$ the predefined semiaxes length of the ellipsoid. The data points $\{p_{i}\},\ i=(1,2,\ldots,N\geq 6)$ on the surface of this ellipsoid are generated most conveniently in polar coordinates:

[TABLE]

where the angles $\theta_{i}$ and $\phi_{i}$ are generated from uniformly distributed random numbers in the interval $[0,\pi]$ and $[0,2\pi]$ respectively. $\bf R$ is the rotation matrix whose elements are given in terms of the Euler angles as:

[TABLE]

We generate several sets of data for aligned as well as arbitrarily oriented ellipsoids corresponding to different ellipsoidal shapes characterised by a parameter $\chi=A/C$; defined as the ratio of major and minor axes. Different data sets were generated corresponding to $\chi$ values ranging from $\chi\sim 1.5$ to $\chi\sim 10^{4}$ to test this algorithm. A particular case of interest are extremely elongated and flat ellipsoids with a very high value of $\chi$ - such ellipsoids arise in the context of gravitational wave data analysis and are separately discussed in the next section.

Results

Some sample results of the ellipsoid reconstruction tests using synthetic data sets are summarised in the tables below.

The first column of each table are input semi-axes length and Euler angles $A,B,C,\alpha,\beta,\gamma$ used to generate the data points using Eq (18). The second and third columns are the estimated values of these quantities using the algorithm described in this paper. As mentioned earlier, there are two independent ways to make the initial guess for $\mathbf{R}$ used to project the data. One involves initial estimation from the eigenvectors of the Fisher information matrix (calculated from input data) whereas the other involves initializing $\mathbf{R}$ using a random positive definite matrix whose elements are drawn from $\mathcal{N}(0,1)$. The tables show the result of both these cases. The fourth column contains the estimate of ellipsoid parameters using Li’s algorithm and serves as a baseline. Data sets 1-3 tabulate sample results where the input data is generated from an ellipsoidal surface whose principal axes are aligned with $XYZ$. Data sets 4-6 are sample results for the non-aligned case.

For aligned cases corresponding to small values of $\chi$ (e.g. Data Set-1), both algorithms (the one presented in this paper and Li’s (2004)) were able to reconstruct the ellipsoid accurately. But for $\chi$ values $\geq 5$ (e.g. Data Sets-2, 3, 5, 6), Li’s method was observed to reconstruct the ellipsoid incorrectly in certain cases. This issue becomes more pronounced as we increase $\chi$. For non-aligned ellipsoids with $\chi\sim 10^{2}$ or more, Li’s method gave different answers for every new set of random input points. On the other hand, the algorithm presented here was robust even for extreme ellipsoidal shapes arbitrarily oriented to the fixed axes.

The algorithm as given by Ying X et al. [9] works robustly for ellipsoidal surfaces whose principal axes are aligned along the fixed frame of reference. It is seen to perform well for extremely elongated and flat ellipsoids (large values of $\chi$). But for arbitrarily oriented ellipsoids, it was observed that the method was not able to reconstruct the ellipsoids correctly. Therefore the method is very sensitive and numerically unstable for non-aligned ellipsoidal surfaces.

V Case study: Application to Gravitational Wave Data Analysis

Gravitational wave data analysis aims to search for weak gravitational wave signals from compact massive astrophysical objects such as binary neutron stars and black holes. The direct detection of astrophysical gravitational wave signals has opened a new observational window to the physical Universe, and is expected to complement the information obtained from traditional astronomy in various electromagnetic observation bands.

The two advanced LIGO [13] gravitational-wave observatories made the first direct detection of these signals about a year ago on 14 September 2015 [14]. Several other observatories are also being commissioned including LIGO-India [15]. A network of three or more of such detectors is expected to improve the scientific potential of these searches by a large margin.

Accurate gravitational wave signal from inspiral, merger and coalescence of compact binaries can be calculated [16] theoretically. These theoretical models allow the well-known technique of matched filtering to be used for detecting faint signals buried in detector noise. If the latter is Gaussian, then it can be shown that the matched filter is optimum, yielding maximum signal to noise ratio (SNR).

The matched filtering strategy for gravitational wave searches is to compute the cross-correlation between the interferometer output and a set of template waveforms; over the detector bandwidth, weighted inversely by the noise power spectrum of the detector ([17, 18, 19]). The deemed parameter space is gridded appropriately for adequate coverage, and the template waveforms are constructed for every point in this grid. The grid is also known as the bank of templates in gravitational wave literature. The construction of this template bank (in other words, the grid over the deemed parameter space) is aided by inducing a metric on the signal manifold.

The match, or overlap between two templates $h(\vec{\lambda}_{1})$ and $h(\vec{\lambda}_{2})$ is defined through their inner-product integral

[TABLE]

where, $\sim$ over the symbols denotes frequency domain representation and $*$ denotes complex conjugation. The one-sided noise power spectral density is given by $S_{h}(f)$.

For two nearby templates $h(\vec{\lambda})$ and $h(\vec{\lambda}+\Delta{\vec{\lambda}})$ one can expand the match $M$ as a Taylor series around $\vec{\lambda}$, to get the following expression upto lowest order term in $\Delta{\vec{\lambda}}$:

[TABLE]

In the above expression, $\Delta\lambda^{i,j}\in\Delta\vec{\lambda}$. Further, we have normalized the templates such that $M(\vec{\lambda},\Delta{\vec{\lambda}}=0)=1$. From Eq. (21) one can see that the $\it{mismatch}$ $(1-M)$ can be used to define the distance square $(\Delta s^{2})$ between two nearby templates in terms of a metric $g_{ij}$ induced on the signal manifold as:

[TABLE]

where,

[TABLE]

It is clear that at a fixed minimal match, the above equation describes the surface of a hyper-ellipsoid (fixed centre).

For a $N$ dimensional signal manifold, $g_{ij}$ is a square symmetric matrix with $N(N+1)/2$ independent components. We pause to note that $g_{ij}$ is not constant over the signal manifold due to to the curvature of the space. This metric is widely used in the gravitational wave signal analysis for placement of templates, determination of consistency of triggers from multiple detectors [20], etc.

Our aim in this section is to demonstrate the efficient numerical estimation of $g_{ij}$ using the technique developed in earlier sections of this paper for the particular case of $N=3$. The latter corresponds to the case where the signal is described by three parameters: two component masses and an effective mass weighted spin magnitude parameter [21] of the compact binary system. Conventionally, one reparameterizes these to new chirp time coordinates $(\tau_{0},\tau_{3},\tau_{3s})$ in which the metric is almost flat (slowly changing). Comparing Eq. (23) with Eq. (13), we immediately notice the correspondence between the metric $g_{ij}$ and the Fisher information matrix $\mathbf{K}$. In this case $g_{ij}$ at a fixed point $\vec{\lambda}$ has six independent components. One starts by numerically solving the match equation $\left\langle h(\vec{\lambda}),\;h(\vec{\lambda}+\Delta{{\vec{\lambda}}})\right\rangle=M$ for $\Delta\vec{\lambda}$ along (at least) 6 random directions. These six points can now be used to fit the ellipsoidal surface Eq. (23) using the iterative technique developed earlier. The Fisher information matrix corresponding to the “best-fit” ellipsoid gives the best numerical estimate of the metric $g_{ij}$ at the point $\vec{\lambda}$ in the parameter space.

We demonstrate the ellipsoid fitting and numerical estimation of the metric $g_{ij}$ on the gravitational wave signal manifold, at a point $\vec{\lambda}_{0}$ correspond to the component masses $(8.0,2.0)M_{\odot}$ and reduced spin magnitude $0.0$. The minimal match is taken to be $0.97$ and the advanced LIGO ’aLIGOZeroDetHighPower’ model for $S_{h}(f)$ is assumed. The match equation is solved along six randomly chosen directions from $\vec{\lambda}_{0}$ and used as test data points to fit the ellipsoid which is shown in Fig 3. Table-VII shows the estimated semi-axes and Euler angles (in degrees) obtained from $g_{ij}$:

We notice that the ellipsoid has high value of $\chi\simeq\frac{2.14}{.004}=535$. The elements of the Fisher information matrix can also be estimated using semi-analytic techniques, using the classic formula

[TABLE]

We find that the numerical approach via ellipsoid fitting provides a better estimate of the metric as determined by the following test. The metric components are calculated at a fiducial point ($80.0,3.5,0.0$) in the chirp time ($\tau_{0},\tau_{3},\tau_{3s}$) coordinate space at a constant minimal match of $M=0.97$ in two different ways: (a) by semi-analytical technique outlined in Eq. (25) and (b) numerically, via the ellipsoid fitting algorithm presented in this work. If the metric coefficients are correctly computed, then we expect the squared distance between a random point on the ellipsoidal surface and the centre to be $\leq(1-M)$. In other words, an arbitrary point on the ellipsoids should have a match $\geq M$ with the centre as evident from Eq. (24). To this end, we sprayed $40,000$ points at random on both these ellipsoidal surfaces and plotted their match with the centre in Fig. 4. The colour bar in the figure indicates the numerical range of matches. It is clear that the metric estimated by ellipsoidal fitting (left panel) is better in comparison to the semi-analytical method, as the matches are closer to the desired level (0.97). On the other hand, many parts of the semi-analytically obtained ellipsoid show match values that are significantly below this level (indicating over-estimation of the ellipsoid size).

VI Conclusion

In this paper, we have developed a general algorithm for fitting ellipsoids of arbitrary shape and orientation using an iterative random transformation based method. We have shown the new method can fit long, thin or compressed ellipsoid and also able to retrieve the rotation angle accurately. Our method is based on iteratively improving the fit by changing the orientation of the coordinates to align along the axes of the ellipsoid. We have verified the accuracy and numerical stability of our algorithm using several sets of synthetic data. Finally, we have demonstrated how this algorithm can be used to numerically estimate the metric on the signal manifold of gravitational wave signals.

Acknowledgements

Amit Reza would like to thank Indian Institute of Technology Gandhinagar for a research fellowship. Thanks are also due to fellow PhD students Soumen Roy, Chakresh Kr. Singh and Md. Yousuf Jamal for useful discussions and help with the manuscript.