Finding the orientation of a noisy ellipse

Question

This question comes from a neuroscience study which generates $12$ vectors. The vectors are evenly spaced, $30 n$ degrees for $n=0,\dots, 11$, each with their tail centered on the origin.

I am looking for biases in each set of vectors, in which one orientation is favored. By this I mean that there may be $2$ peaks opposed from one another (by $180^{\circ}$). So I am wondering what you think is the simplest way to estimate the "orientation" with the greatest magnitude, given that the measurements are discrete, and the true peak may be between datapoints.

The data can be considered random, but the ideal case would look something like an ellipse centered on the origin, with $2$ equal and opposed maxima and minima. I suppose another way to frame the question would be how to find the direction of the major axis of an ellipse when plotted discretely over $12$ uniformly distributed theta values between $0^{\circ}$ and $360^{\circ}$. But one has to consider that the plot may also include significant random noise, so finding one peak is not enough to definitively tell the orientation of the ellipse.

I first learned to do this for a single peak, which would resemble a cardioid. This was simple: the direction of bias could be estimated by taking the vector sum of the $12$ vectors. Now that I am working with $2$ peaks, I do not feel as certain about my technique, but I have made an attempt:

• Take the sum of the magnitudes of opposing vectors

• Assign them to 6 evenly spaced vectors with $\theta = {60n}^{\circ}$ for $n=0,\dots, 11$

• Vector sum

• Divide $\theta$ of $V_{sum}$ by $2$

This solution seems too simple to be right... and I don't have the skill to prove or disprove it. Would this reliably determine the orientation of a low resolution ellipse?

Answer 1

Consider the coordinates as complex numbers; let u+iv be the sum of their squares; take ½ arctan(v/u).

Later: Or rather, atan2(v,u)/2; if u<0, the naïve form above will give you the short axis instead. In Python, I'd use cmath.phase(z)/2.

Much later: I should have said in the first place that this is how to find the line through the origin that minimizes the squared distance of the sample points to the line. If the samples lie on an ellipse centred at the origin, it stands to reason that this line is close to the ellipse's long axis. But if the axis is not one of the sample lines, there may be a bias; one of these years I'll look into that.

Answer 2

I would build an outer product tensor and find it's eigensystem. First to build the tensor: $${\bf T} = \sum_{samples}\left(\sum_{i=0}^{11} {\bf v}_i {\bf v}_i^T\right)$$

By the spectral theorem, since it will be symmetric (why?), $\bf T$ is then ensured to have ON system with eigenvalues $\lambda_k > 0$ and normalized eigenvectors ${\bf \hat e}_k$ $${\bf T} = \sum_{i=1}^{\dim({\bf T})} \lambda_i {\bf \hat e}_i{\bf \hat e}_i^T$$ Since they are real (non-negative even), we can sort them and any orientational bias would be if they differ from each other. You can measure this in many ways, maybe variance of eigenvalues, or just difference between smallest and largest and so on.

---

The benefit of using this method rather than many of the others proposed here is that it easily expands to any number of dimensions. Fourier transforms on spheres are tedious, complex numbers can treat two dimensions but already at three it become a nuisance and so on.

Answer 3

If the experimental points appear roughly on an ellipse (of course with scatter), you can use a regression technic to compute the parameters of the theoretical ellipse. With the equation obtained, it is easy to determine the axes of the ellipse and check how far is the center of the ellipse from the origine.

This comes from the paper : https://fr.scribd.com/doc/14819165/Regressions-coniques-quadriques-circulaire-spherique , page 16.

For more information about the parameters and properties of the ellipses: http://mathworld.wolfram.com/Ellipse.html

If we want to fit an ellipse centred at the origine the matrix system is reduced to :

Answer 4

I would assume that your points $(x_i,y_i)$, $i=1,\dots, N$ come from some normal bivariate distribution with mean $\bar x= 0$ and $\bar y = 0$. In such case you can easily find the covariance matrix: $$ \Sigma = \begin{pmatrix}a & b \\ b & c\end{pmatrix} $$ with these formulas: $$ a = \frac 1 N \sum_{i=1}^N x_i^2, \qquad b = \frac 1 N \sum_{i=1}^N x_i y_i, \qquad c = \frac 1 N \sum_{i=1}^N y_i^2. $$ Then you should compute the two eigenvectors of $\Sigma$ (or $\Sigma^{-1}$? ...to be checked) and corresponding eigenvalues. The direction you look for should be the eigenvector corresponding to the largest eigenvalue.

(see https://stats.stackexchange.com/questions/24380/how-to-get-ellipse-region-from-bivariate-normal-distributed-data)

Answer 5

Take the Discrete Fourier Transform of the vector lengths $\bf v$.

The 2nd order coefficient, actually its phase, gives you the angle you are requesting; remember to divide by two.

Octave (Matlab) code is straightforward:

F=fft(V);
a=angle(F(3))/2;   # 3rd element = 2nd order coefficient

but with 12 points you can easily write down the full expansion for the desired coefficients as a linear combination of all vectors where coefficients are sines and cosines of easy angles (and you can obtain the real/imag parts separately with operations that just involve real numbers).

Then you only need one transcendent operation i.e. the arctan() to get the angle.

The strength of the method is, that the result is invariant for a shift of the 12 vector pattern, as required.

---

But there's more. Beside the requested angle, it gives you a measurement of the mentioned "noise".

In fact, the ratio $$ \frac{2 \cdot |F(3)|^2}{12 \cdot \sum{\bf v}_i^2} (*) $$ measures how much the two-lobe effect (this is what is meant with "ellipse", afaik) is relevant compared to the rest of the noise. You can use this ratio as a percentage, as it is energetically meaningful (by Parseval's theorem)

(*) Note. The $12$ factor at denominator depend on actual formula used for discrete Fourier transform. The above is valid for Octave (Matlab).

More in detail, using also the other coefficients, you can tell how big is the "two-lobe" effect compared to various k-lobed "modes" (you would probably call noise only those with k>2).

In fact, once found the F coefficients, their modula measure the amplitudes of the n-lobe modes. I.e. more exactly:

• 0th order i.e. $|F(1)|^2$, gives the weight of the constant mode (0-lobes)
• kth order, $|F(1+k)|^2+|F(12-k)|^2 == 2 \cdot|F(1+k)|^2$, is the weight of the k-lobe mode (k in 1:5)
• 6th order, $|F(7)|^2$, gives the weight of the the 6-lobed mode, i.e. the maximum angular frequency mode.

In conclusion:

• $\frac{2 \cdot |F(3)|^2}{12 \cdot \sum{\bf v}_i^2}$ is the percentage of the data energy which is reated to the two lobes
• $\frac{|F(1)|^2 + 2 \cdot |F(2)|^2}{12 \cdot \sum{\bf v}_i^2}$ measures the low frequency "noise" (or bias, name it as you want)
• $\frac{2 \cdot |F(4)|^2+2 \cdot |F(5)|^2+2 \cdot |F(6)|^2 + |F(7)|^2}{12 \cdot \sum{\bf v}_i^2}$ measures the energy in the higher order modes (call it the hi frequency noise)

That's why I find the Fourier representation of variability particularly meaningful when explaining variability of data representable in polar plots.