# Finding the orientation of a noisy ellipse

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?

• Keep $\theta=30°n$. – Yves Daoust Sep 14 '15 at 8:04
• @Ryan : I not fully see how your graph looks like.Can you post an example. – JJacquelin Sep 14 '15 at 8:05
• Are your data points truly aligned on an ellipse (with added noise), or just on some elongated curve ? – Yves Daoust Sep 14 '15 at 8:11
• The bottom right of this image is the ideal scenario, so it is not a true ellipse. However, the rest of the plots in this image are also common outcomes, and it is expected that there will be a range. I just need to be able to find the orientation, in the event that it resembles an ellipse. – Ryan Sep 14 '15 at 8:19
• These plots are not at all ellipses. You should discard this term from the problem statement, as it is completely misleading. From what I see, in several cases associating a dominant direction to a plot is a nonsense. – Yves Daoust Sep 14 '15 at 8:57

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.

• Thank you, I have already implemented this method and it works well. I have to accept your response, as it is certainly the simplest, and it does exactly what I need it to. If you have the time, I would really appreciate your help in understanding how you arrived at your solution. Here are some of the graphs I have generated using your method (in yellow), mine (in blue), and the standard unidirectional vector sum (red) – Ryan Sep 14 '15 at 23:24
• Awhile ago I needed a least-squares fit of a line to a given set of noisy points: not treating y as a linear function of x, but minimizing the orthogonal distance to the line. I did the algebra and this is what I got (after shifting the origin to the average). I treated your question as equivalent to this: if the points do lie on an ellipse, obviously the long axis minimizes their distances. – Anton Sherwood Sep 14 '15 at 23:35
• Informally: notice that squaring a complex number brings the images of $+z$ and $-z$ together; and, of course, the longest vectors contribute most to the sum. Thus the sum of complex squares is the square of a number representing your long axis. The negation of that sum is the square of a number representing the short axis, a right angle away. – Anton Sherwood Sep 14 '15 at 23:39
• Modifying my first comment: “if the points do lie on an ellipse” should be “if the points are symmetrically distributed on an ellipse”. I wouldn't have proposed this rule if you hadn't specified that the given points are nicely distributed. – Anton Sherwood Sep 14 '15 at 23:43
• Almost forgot: Thanks for sharing your results! – Anton Sherwood Sep 14 '15 at 23:46

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.

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.

• Since we only need one component fft is a bit overkill. Just multiplying with the third complex exponential frequency is $\mathcal{O}(n)$ slightly better than fft $\mathcal{O}(n\log(n))$ for calculating the whole fft. Still a good answer. Don't know why you got a down vote for it. – mathreadler Sep 15 '15 at 2:11
• Oh, wait. I see you mentioned that later. Nevermind. – mathreadler Sep 15 '15 at 2:17
• The FFT offers many more advantages (to me at least). tried to insert some hints in the post. And with n=12, the O(n) issue is really not an issue, I believe.I don't understand the reason of the down vote either. But it's my first answer. Where can I see motivations, if any? Anyway, I'll improve. – lurix66 Sep 15 '15 at 12:11
• People don't need to give reasons for down-votes which can be frustrating sometimes as it's difficult to know how to improve without feedback. It is recommended to learn Latex for typesetting. I can show you how by helping to edit this answer. – mathreadler Sep 15 '15 at 12:25
• Now if you accept my edit and hit "edit" you can see how you can embed latex math formatting syntax into answers and questions. Sometimes people get frustrated and downvote just because stuff is unreadable. – mathreadler Sep 15 '15 at 12:35

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.

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.