8

I have the following data:-

  • I have two points ($P_1$, $P_2$) that lie somewhere on the ellipse's circumference.
  • I know the angle ($\alpha$) that the major-axis subtends on x-axis.
  • I have both the radii ($a$ and $b$) of the ellipse.

I now need to find the center of this ellipse. It is known that we can get two possible ellipses using the above data.

I have tried solving this myself but the equation becomes so complex that I always give up.

This is what I have done till now:-

I took the normal ellipse equation $x^2/a^2 + y^2/b^2 = 1$. To compensate for the rotation and translation, I replaced $x$ and $y$ by $x\cos\alpha+y\sin\alpha-h$ and $-x\sin\alpha+y\cos\alpha-k$, respectively. $h$ and $k$ are x and y location of the ellipse's center.

Using these information I ended up with the following eq:- $$a B_1\pm\sqrt{a^2 B_1^2 - C_1(b^2 h^2 - 2 A_1 b^2 h)} = a B_2\pm\sqrt{a^2 B_2^2 - C_2 (b^2 h^2 - 2 A_2 b^2 h)} \quad (1)$$

where $A = x\cos\alpha +y\sin\alpha$, $B = -x\sin\alpha+y\cos\alpha$ and $C = a^2 B^2 + A^2 b^2 - a^2 b^2$.

Now the only thing I need to get is $h$ from (1). All other values are known, but I am not able to single that out.

Anyway if the above equations looks insane then please solve it yourself, your way. I could have drifted into some very complicated path.

AppleGrew
  • 257
  • @Peter: But the subsequent equations 21 and 22 make it seem difficult to integrate the given information about the semi-axes into that approach. – joriki Jul 22 '11 at 14:21
  • I was hoping this would involve using an eccentric anomaly but seems you guys have solved without :) – Alice Jul 22 '11 at 16:59

4 Answers4

9

Let the points be $P_1(x_1, y_1)$ and $P_2(x_2, y_2)$, assumed to lie on an ellipse of semiaxes $a$ and $b$ with the $a$ axis making angle $\alpha$ to the $x$ axis.

Ellipse

Following @joriki, we rotate the points $P_i$ by $-\alpha$ into points

$$Q_i(x_i \cos(\alpha) + y_i \sin(\alpha), y_i \cos(\alpha) - x_i \sin(\alpha)).$$

We then rescale them by $(1/a, 1/b)$ to the points

$$R_i(\frac{x_i \cos(\alpha) + y_i \sin(\alpha)}{a}, \frac{y_i \cos(\alpha) - x_i \sin(\alpha)}{b}).$$

Circle

These operations convert the ellipse into a unit circle and the points form a chord of that circle. Let us now translate the midpoint of the chord to the origin: this is done by subtracting $(R_1 + R_2)/2$ (shown as $M$ in the figure) from each of $R_i$, giving points

$$S_1 = (R_1 - R_2)/2, \quad S_2 = (R_2 - R_1)/2 = -S_1$$

each of length $c$. Half the length of that chord is

$$c = ||(R_1 - R_2)||/2 = ||S_1|| = ||S_2||,$$

which by assumption lies between $0$ and $1$ inclusive. Set

$$s = \sqrt{1-c^2}.$$

The origin of the circle is found by rotating either of the $S_i$ by 90 degrees (in either direction) and rescaling by $s/c$, giving up to two valid solutions $O_1$ and $O_2$. (Rotation of a point $(u,v)$ by 90 degrees sends it either to $(-v,u)$ or $(v,-u)$.) For example, in the preceding figure it is evident that rotation $R_1$ by -90 degrees around $M$ and scaling it by $s/c$ will make it coincide with the circle's center. Reflecting the center about $M$ (which gives $2M$) produces the other possible solution.

Unwinding all this requires us to do the following to the $O_i$:

  • Translate by $(R_1+R_2)/2$,
  • Scale by $(a,b)$, and
  • Rotate by $\alpha$.

The cases $c \gt 1$, $c = 1$, and $c=0$ have to be treated specially. The first gives no solution, the second a unique solution, and the third infinitely many.

FWIW, here's a Mathematica 7 function. The arguments p1 and p2 are length-2 lists of numbers (i.e., point coordinates) and the other arguments are numbers. It returns a list of the possible centers (or Null if there are infinitely many).

f[\[Alpha]_, a_, b_, p1_, p2_] := Module[
   {
    r, s, q1, q2, m, t, \[Gamma], u, r1, r2, x, v
    },
   (* Rotate to align the major axis with the x-axis. *)
   r = RotationTransform[-\[Alpha]];
   (* Rescale the ellipse to a unit circle. *)
   s = ScalingTransform[{1/a, 1/b}];
   {q1, q2} = s[r[#]] & /@ {p1, p2};
   (* Compute the half-length of the chord. *)
   \[Gamma] = Norm[q2 - q1]/2;
   (* Take care of special cases. *)
   If[\[Gamma] > 1, Return[{}]];
   If[\[Gamma] == 0, Return[Null]];
   If[\[Gamma] == 1, 
    Return[{InverseFunction[Composition[s, r]][(q1 + q2)/2]}]];
   (* Place the origin between the two points. *)
   t = TranslationTransform[-(q1 + q2)/2];
   (* This ends the transformations.  
   The next steps find the centers. *)
   (* Rotate the points 90 degrees. *)
   u = RotationTransform [\[Pi]/2];
   (* Rescale to obtain the possible centers. *)
   v = ScalingTransform[{1, 1} Sqrt[1 - \[Gamma]^2]/\[Gamma]];
   x = v[u[t[#]]] & /@ {q1, q2};
   (* Back-transform the solutions. *)
   InverseFunction[Composition[t, s, r]] /@ x
   ];
whuber
  • 8,531
  • Wow! Thanks. Will try and understand this. +1 for now. – AppleGrew Jul 22 '11 at 20:21
  • Here's the x-coordinate ;-) $\frac{1}{2 a b}\left(a b (x_1+x_2)+\left(a^2 (-y_1+y_2) \cos[\alpha]^2+(a-b) (a+b) (x_1-x_2) \cos[\alpha] \sin[\alpha]+b^2 (-y_1+y_2) \sin[\alpha]^2\right) \surd \left(\left(b^2 \left((x_1-x_2)^2+(y_1-y_2)^2\right)+a^2 \left(-8 b^2+(x_1-x_2)^2+(y_1-y_2)^2\right)+(a-b) (a+b) (-(x_1-x_2+y_1-y_2) (x_1-x_2-y_1+y_2) \cos[2 \alpha]-2 (x_1-x_2) (y_1-y_2) \sin[2 \alpha])\right)/\left(-\left(a^2+b^2\right) \left((x_1-x_2)^2+(y_1-y_2)^2\right)+(a-b) (a+b) ((x_1-x_2+y_1-y_2) (x_1-x_2-y_1+y_2) \cos[2 \alpha]+2 (x_1-x_2) (y_1-y_2) \sin[2 \alpha])\right)\right)\right)$ – whuber Jul 22 '11 at 21:54
  • are you sure this equation is correct? I am not getting the correct results. And I am guessing $\sin[\alpha]^2$ means sine of square of alpha. – AppleGrew Jul 25 '11 at 02:40
  • You should check http://jsfiddle.net/6Wxbg/ . The two pink points are the given points. Ideally the two ellipse must intersect at the two points. – AppleGrew Jul 25 '11 at 04:18
  • @J.M. "Here's" where? – AppleGrew Jul 25 '11 at 04:43
  • Sorry, I hit the Enter key too early; let me get back to you... – J. M. ain't a mathematician Jul 25 '11 at 04:45
  • @AppleGrew: Strange, you write "And I am guessing $\sin[\alpha]^2$ means sine of square of alpha", but your code correctly calculates the square of the sine of $\alpha$. – joriki Jul 25 '11 at 05:13
  • @Joriki I modified it after J.M.'s response. His formula clearly states that it should be square of the sin and cos. – AppleGrew Jul 25 '11 at 09:11
  • @whuber: I edited your math display because what you typed is too long and is breaking the webpage. Please consider editing your expression in to your answer in a sensible way taking care of line breaks. – Willie Wong Jul 25 '11 at 13:31
  • @Willie Thanks much. The point of that comment lies in the length and complexity of the formula, which clearly is inferior (both conceptually and computationally) to the algorithm presented in the reply itself. In short, I view the line overflow as part of the message, not a formatting error. – whuber Jul 25 '11 at 13:48
  • @whuber I never tried to understand your solution. I was trying to eat up cooked meal, but I now understand, I will have to cook my own JS algo using your ingredients. I now have some doubts. I guess $(R_1+R_2)/2$ means $x=(x_1+x_2)/2$ and $y=(y_1+y_2)/2$. Then I don't understand, how $S_2=-S_1$? From your solution I get the feel that you have assumed that chord is parallel to the axes. Is that true? If yes then that's a problem as that will not be the case here. We 'un-rotate' only the ellipse not the chord itself. – AppleGrew Jul 25 '11 at 18:34
  • @Apple Nope, no assumptions: the solution is completely general and covers all possibilities. Yes, you're correct: $(R_1+R_2)/2$ means average the coordinates. Once you move the origin there, the two points are necessarily negatives of each other. – whuber Jul 25 '11 at 20:19
  • @Whuber If possible can then can you provide a diagram? I guess I have the wrong diagram in my head. – AppleGrew Jul 25 '11 at 20:49
  • @AppleGrew Done. – whuber Jul 25 '11 at 21:28
  • @whuber Superb!But if I understand it correctly then the fig. 2 has some problem. This is a 4 stage process. Stage1 scale to unit circle. Stage2 translate chord mid-pt to origin. Stage3 rotate and scale to get the origin. Stage 4 invert all. If fig.2 shows the state after stage1 then the circle shouldn't necessarily be at origin. If it shows after stage2 then $M$ should be at origin and chord ends should be $S_1$ and $S_2$. Also isn't expressing $c$ as $c=||S_i||/2$ more appropriate? I got confused by $||S_1||$, as that denotes distance of point $S_1$ from origin, i.e. half of chord length. – AppleGrew Jul 26 '11 at 06:00
  • @Apple It doesn't matter where the circle is situated. Once the problem has been transformed into the question of finding a unit circle containing two points, it has many solution methods. Mine simply rotates the points 90 degrees about their midpoint and rescales them. There exist many, many other solutions, beginning with simple compass and ruler constructions of the ancient Greeks. The points $S_i$, by the way, are points in a coordinate system in which $M$ has been moved to the origin, so indeed $||S_i||$ = $c$ = $||R_2-R_1||/2$, as shown. – whuber Jul 26 '11 at 13:35
  • [Updated] To All. A demo of this code implemented using Javascript can be found at http://jsfiddle.net/ZxRBT/. – AppleGrew Sep 19 '11 at 06:26
  • I tried this (and the example implemented in JS @AppleGrew) and already the first step does not work. When you rotate the points and scale them back, the resulting points are > 1 (no union circle). If you change any parameter in the JS example the center is not drawn correctly anymore – Fuzzyma Oct 24 '15 at 22:03
  • @whuber This new question may also involve a formula for the center of the ellipse. Can you kindly take a look? – Tito Piezas III Jan 15 '24 at 14:51
  • @whuber I have a another question about pairs of lattice point triplets on an ellipse. Might be of interest. – Tito Piezas III Jan 16 '24 at 08:33
  • @Tito Thanks, but that looks like it's really a question about rational points on elliptic curves. – whuber Jan 16 '24 at 15:04
  • @whuber Elliptic curves may be involved, but not totally. We require the equation $x^3+y^3+z^3 = (z+1)^3$ to be solved in the integers, not in the rationals. – Tito Piezas III Jan 16 '24 at 15:14
  • @Tito It amounts to the same thing. One usually projectivizes the equation, resulting in rational solutions, and finds the integral solutions among them. In doing so, you can exploit the group of the curve. – whuber Jan 16 '24 at 15:17
  • 1
    @whuber Ok. I guess I better add the tag "elliptic curves" then. – Tito Piezas III Jan 16 '24 at 15:18
7

You can make your life a lot easier by rotating the two given points through $-\alpha$ instead of rotating the coordinate system through $\alpha$. Then you can work with the much simpler equation

$$\frac{(x-x_0)^2}{a^2}+\frac{(y-y_0)^2}{b^2}=1\;.$$

If you substitute your two (rotated) points $(x_1,y_1)$ and $(x_2,y_2)$, you get two equations for the two unknowns $x_0$,$y_0$. Subtracting these from each other eliminates the terms quadratic in the unknowns and yields the linear relationship

$$\frac{x_1^2-x_2^2-2(x_1-x_2)x_0}{a^2}+\frac{y_1^2-y_2^2-2(y_1-y_2)y_0}{b^2}=0\;.$$

You can solve this for one of the unknowns and substitute the result into one of the two quadratic equations, which then becomes a quadratic equation in the other unknown that you can solve.

Of course in the end you have to rotate the centre you find back to the original coordinate system.

joriki
  • 238,052
  • 5
    +1 Even easier: after the rotation, rescale the axes by $1/a$ and $1/b$. This reduces the problem to finding the center of a unit circle given two points on it, which is a standard (easy) Euclidean geometry construction. – whuber Jul 22 '11 at 15:10
  • Ah, yes, that makes sense :-) – joriki Jul 22 '11 at 15:12
  • 3
    Indeed, you can then further rotate and translate the configuration to place the two points at $(0,\pm c)$. Clearly then the solutions for the center are $(\pm s, 0)$ where $s^2=1-c^2$. (Obviously there's no solution when $|c| \gt 1$.) – whuber Jul 22 '11 at 15:21
  • You should write that as an answer; it's better than mine :-) – joriki Jul 22 '11 at 15:26
  • It's the same as yours; these are just comments :-). – whuber Jul 22 '11 at 15:58
  • @whuber: The affine transformation trick for ellipses is really the gift that keeps on giving... :) – J. M. ain't a mathematician Jul 22 '11 at 16:07
  • Sure: all ellipses are just the unit circle. If the problem also permits the application of projective transformations, all (nondegenerate) conic sections are also just the unit circle. That immediately lets us solve a generalization of this problem to hyperbolae and parabolae. – whuber Jul 22 '11 at 16:48
  • Ooh! This looks promising. Will try to evaluate this once more. – AppleGrew Jul 22 '11 at 17:29
  • Aah! One question. How do I find the the converted P1 and P2? – AppleGrew Jul 22 '11 at 18:18
  • You rotate them just like you rotated $x$ and $y$, just with the opposite angle. – joriki Jul 22 '11 at 18:19
  • @whuber would you spoon feed your answer? I have been trying to solve this since one week. Well the original problem was bigger but then it 'simplified' to this. – AppleGrew Jul 22 '11 at 18:33
  • @joriki The question about the center of the ellipse has been answered. This new question about pairs of lattice point triplets on an ellipse may be of interest though. – Tito Piezas III Jan 16 '24 at 08:31
2

Since the expression in whuber's comment was too darn long, here's the x-coordinate expression:

$$\begin{align*} &\frac1{2ab}\left(ab (x_1+x_2)+\left(a^2 (y_2-y_1)\cos^2\alpha+(a-b)(a+b) (x_1-x_2) \cos\,\alpha\sin\,\alpha+b^2 (y_2-y_1)\sin^2\alpha\right)\right.\\ &\left.\surd \left(\left(b^2 \left((x_1-x_2)^2+(y_1-y_2)^2\right)+a^2 \left(-8 b^2+(x_1-x_2)^2+(y_1-y_2)^2\right)+\right.\right.\right.\\ &\left.\left.\left.(a-b)(a+b) (-(x_1-x_2+y_1-y_2) (x_1-x_2-y_1+y_2) \cos\,2\alpha-2 (x_1-x_2) (y_1-y_2) \sin\,2\alpha)\right)/\right.\right.\\ &\left.\left.\left(-\left(a^2+b^2\right) \left((x_1-x_2)^2+(y_1-y_2)^2\right)+(a-b) (a+b) ((x_1-x_2+y_1-y_2) (x_1-x_2-y_1+y_2) \cos\,2\alpha+\right.\right.\right.\\ &\left.\left.\left.2(x_1-x_2)(y_1-y_2)\sin\,2\alpha)\right)\right)\right) \end{align*}$$

  • Well checkout http://jsfiddle.net/6Wxbg/1/ . As I mentioned this is not giving me the desired result. The two ellipses should intersect at the pink dots. – AppleGrew Jul 25 '11 at 05:04
  • I don't understand why there is "the x-coordinate" -- shouldn't there be two x-coordinates, one for each centre? The green ellipse seems to be on target (it hits the top left corners of the pink dots, which is correct since you draw those from there), so this x-coordinate seems to be the correct one for cy2, and you probably need to change a sign somewhere to make it work for cy1 -- you can sort of tell that the red ellipse will be OK if you slide it to the right. – joriki Jul 25 '11 at 05:30
  • I haven't checked whuber's solution in detail (yet) so I can't say anything... – J. M. ain't a mathematician Jul 25 '11 at 05:32
  • By symmetry, I think it must be the $+$ after $ab(x_1+x_2)$ that needs to be changed to a $-$ for cy1, since the solutions should be symmetric around $(x_1+x_2)/2$. However, there seems to be an additional problem -- when I change y1 to $50$, both ellipses are off. – joriki Jul 25 '11 at 05:41
  • I think there must be something wrong with B -- it's the average of the two y-coordinates of the centres, but it's being calculated only in terms of x0 and y0; that can't be right. – joriki Jul 25 '11 at 06:11
  • @joriki Are you referring to eq $B=y0\cos(ang)-x0\sin(ang)$? Here the constants are not symmetric. I derived this by substituting the rotated and then translated values of x and y and substituted then into standard ellipse equation, as per your suggestion. Then I solved this quadratic equation for $y_c$ (y-coord of center). – AppleGrew Jul 25 '11 at 08:42
  • @AppleGrew: Yes, I was referring to that. I didn't look at the derivation in detail, but just from symmetry, it can't be right that the y-coordinates of the two centres are $B+E$ and $B-E$, so their average is $B$, and $B$ is derived only from the coordinates of one of the points -- the average of the two centres must be the average of the two points, by symmetry. – joriki Jul 25 '11 at 08:58
  • Hmm not sure about this but as far as I remember it was fine. I don't have the derivation right now. Will post it when I am back home. – AppleGrew Jul 25 '11 at 09:08
  • It's fine as long as the two y coordinates are $0$ -- you're not actually testing any of the terms proportional to one of the $y$ coordinates with the example you're using, since they're both $0$ -- you can only see this problem when you make them different. – joriki Jul 25 '11 at 09:13
  • @joriki regarding your concern that there should be more x-coords, I guess, we will get the second x by putting a minus before the square root expression. – AppleGrew Jul 25 '11 at 10:38
  • 3
    @All I could provide the y-coordinate of this solution, along with (x,y) for the second solution, but that's not the point. Isn't it obvious that implementing such lengthy convoluted formulas is inferior to implementing the algorithm as described? But if anybody wants the full solution, just ask Mathematica to carry out a FullSimplify of the generic formula FullSimplify[f[\[Alpha], a, b, {x1, y1}, {x2, y2}] (applying suitable domain restrictions and after removing the If statements to check for special cases.) – whuber Jul 25 '11 at 13:51
  • 1
    @J.M. Thanks for introducing me to the align* construct--it's just what is needed. Incidentally, I tested the Mathematica code (for the function f which produced this formula) by creating a lot of random configurations, applying the code, and checking that it returned the correct center every time. Although something may have happened in the process of copying the formula from Mathematica and reformatting it here--which further illustrates the dangers of implementing long complex formulas in code--, I can vouch for the correctness of the overall approach. – whuber Jul 25 '11 at 21:59
  • @whuber: FWIW, that's where I'd use "common expression" isolation. You know, reduced level of abstraction and all that... ;) – J. M. ain't a mathematician Jul 27 '11 at 02:20
1

Unfortunately, my correction on the answer of whuber has been rejected, but I'll try to explain it here.

The chosen answer is very good, but there is an error. $c$ is supposed to be the distance from point $M$ to point $R_1$. The points $R_1$ and $R_2$ are then converted into $S_1$ and $S_2$, respectively, whereby $S_2 = -S_1$. Each $S_1$ and $S_2$ are of length c, as is correctly stated. So the value of c is actually $$c = ||S_1 - S_2||/2 = ||2S_1||/2 = ||S_1||$$ apposed to $c \neq ||S_1||/2$.

It occurred to me while trying to implement the proposed solution. Using this correction, the formula works as expected.

  • Thank you for pointing that out. My exposition indeed was off by a factor of $1/2$--but the code was correct, and no other part of the solution was affected. – whuber Jun 10 '16 at 13:05