8
$\begingroup$

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.

  • 0
    @Peter: But the subsequent equations 21 and 22 make it seem difficult to integrate the given information about the semi-axes into that approach.2011-07-22
  • 0
    I was hoping this would involve using an eccentric anomaly but seems you guys have solved without :)2011-07-22

4 Answers 4

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    ]; 
  • 0
    Wow! Thanks. Will try and understand this. +1 for now.2011-07-22
  • 0
    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)$2011-07-22
  • 0
    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.2011-07-25
  • 0
    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.2011-07-25
  • 0
    @J.M. "Here's" where?2011-07-25
  • 0
    Sorry, I hit the Enter key too early; let me get back to you...2011-07-25
  • 0
    @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$.2011-07-25
  • 0
    @Joriki I modified it after J.M.'s response. His formula clearly states that it should be square of the sin and cos.2011-07-25
  • 0
    @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**.2011-07-25
  • 0
    @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.2011-07-25
  • 0
    @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.2011-07-25
  • 0
    @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.2011-07-25
  • 0
    @Whuber If possible can then can you provide a diagram? I guess I have the wrong diagram in my head.2011-07-25
  • 0
    @AppleGrew Done.2011-07-25
  • 0
    @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.2011-07-26
  • 0
    @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.2011-07-26
  • 0
    [Updated] To All. A demo of this code implemented using Javascript can be found at http://jsfiddle.net/ZxRBT/.2011-09-19
  • 0
    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 anymore2015-10-24
6

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.

  • 4
    +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.2011-07-22
  • 0
    Ah, yes, that makes sense :-)2011-07-22
  • 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$.)2011-07-22
  • 0
    You should write that as an answer; it's better than mine :-)2011-07-22
  • 0
    It's the same as yours; these are just comments :-).2011-07-22
  • 0
    @whuber: The affine transformation trick for ellipses is really the gift that keeps on giving... :)2011-07-22
  • 0
    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.2011-07-22
  • 0
    Ooh! This looks promising. Will try to evaluate this once more.2011-07-22
  • 0
    Aah! One question. How do I find the the converted P1 and P2?2011-07-22
  • 0
    You rotate them just like you rotated $x$ and $y$, just with the opposite angle.2011-07-22
  • 0
    @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.2011-07-22
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*}$$

  • 0
    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.2011-07-25
  • 0
    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.2011-07-25
  • 0
    I haven't checked whuber's solution in detail (yet) so I can't say anything...2011-07-25
  • 0
    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.2011-07-25
  • 0
    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.2011-07-25
  • 0
    @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).2011-07-25
  • 0
    @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.2011-07-25
  • 0
    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.2011-07-25
  • 0
    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.2011-07-25
  • 0
    @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.2011-07-25
  • 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.)2011-07-25
  • 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.2011-07-25
  • 0
    @whuber: FWIW, that's where I'd use "common expression" isolation. You know, reduced level of abstraction and all that... ;)2011-07-27
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.

  • 0
    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.2016-06-10